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Many systems in nature and the synthetic world involve ordered arrangements of units on 
two-dimensional surfaces. We review here the fundamental role payed by both the topology 
of the underlying surface and its Gaussian curvature. Topology dictates certain broad fea- 
tures of the defect structure of the ground state but curvature-driven energetics controls the 
detailed structured of ordered phases. Among the surprises are the appearance in the ground 
state of structures that would normally be thermal excitations and thus prohibited at zero 
temperature. Examples include excess dislocations in the form of grain boundary scars for 
spherical crystals above a minimal system size, dislocation unbinding for toroidal hexatics, 
interstitial fractionalization in spherical crystals and the appearance of well-separated discli- 
nations for toroidal crystals. Much of the analysis leads to universal predictions that do not 
depend on the details of the microscopic interactions that lead to order in the first place. 
These predictions are subject to test by the many experimental soft and hard matter systems 
that lead to curved ordered structures such as colloidal particles self-assembling on droplets 
of one liquid in a second liquid. The defects themselves may be functionalized to create lig- 
ands with directional bonding. Thus nano to meso scale superatoms may be designed with 
specific valency for use in building supermolecules and novel bulk materials. Parameters such 
as particle number, geometrical aspect ratios and anisotropy of elastic moduli permit the 
tuning of the precise architecture of the superatoms and associated supermolecules. Thus the 
field has tremendous potential from both a fundamental and materials science/supramolecular 
chemistry viewpoint. 
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1. Order and geometry in condensed matter 
1.1. Introduction 

More than 200 years ago, in his treatise on the resistance of fluids, d'Alembert 
wrote: "Geometry, which should always follow physics when used to describe na- 
ture, sometimes commands it" [1 . Since then the physics community has explored 
the power of geometry not only to describe, but also to explain structures and 
their properties. In the past 20 years soft condensed matter physics has provided 
many examples of how the geometry of matter is not a quiescent background for 
some microscopic degrees of freedom, but instead plays a major role in determining 
structural and mechanical properties and designing the phase diagram of materials 
such as colloids, liquid crystals, membranes, glasses and carbon nanostructure. 

Geometric models of condensed matter systems have been developed for a wide 
class of materials since the pioneering work of Bernal and Finney [21 |3l SI El [6] . 
In a series of classic papers they suggested that several properties of liquids have 
their geometric counterpart in randomly packed arrays of spheres. The difference 
in density between the solid and the liquid phase of a simple monoatomic sub- 
stance, for instance, is approximately the same between periodically and randomly 
packed hard spheres (roughly 15% — 16%). Also the radial distribution function 
of randomly packed spheres corresponds well with that determined by X-ray and 
neutron diffraction for rare-gas liquids. After Bernal, a significant amount of work 
has been done on random close packing and, even though the legitimacy of the no- 
tion of random close packing itself has been questioned frequently in recent years 
[7], it is now established that many features of the liquid state have in fact a purely 
geometrical nature. 

After the discovery of icosahedral order in metallic glasses jH [9l [10] the idea of 
geometrical frustration (the geometric impossibility of establishing a preferred local 
order everywhere in space, see Sec.|2]), became a fundamental concept for the char- 
acterization of amorphous solids. Farges and coworkers |1H [T2l [T3] were the first 
to show, by electron diffraction experiments and computer simulations, that the 
first atoms of small aggregates of rare gases condensed in ultra-high vacum form 
regular tetrahedra, which later organize in the form of small icosahedral clusters. 
Since icosahedra do not fill three-dimensional Euclidean space R^, the structure 
resulting from the aggregation of these icosahedral building blocks does not ex- 
hibit long range translational order. The lack of crystallization in covalent glasses 
is also rooted in the geometrical frustration associated with the constant coordina- 
tion number of their constituents. Tetravalent monoatomic materials, for instance, 
cannot form a constant angle between bonds incident at the same atom and orga- 
nize in a regular network at the same time. In multiatomic glass-forming materials 
the situation is more involved and the route to the formation of amorphous struc- 
tures is related to the fact that crystallization would require complex activated 
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phenomena and too large a decrease in entropy. 

A breakthrough in the geometrical description of amorphous solids came in 1979 
when Kleman and Sadoc first observed that a number of continuous random lattices 
can be classified as specific mappings of ordered lattices in spaces of constant cur- 
vature onto [14 . This idea was inspired by the observation that whereas regular 
tetrahedra do not fill three-dimensional Euclidean space, they do regularly tile the 
three-dimensional sphere (the manifold described by the equation -^t ~ ^ 
in R^) on which they build a regular polytope (Schlafli symbol {3, 3, 5}) with 120 
vertices of coordination number 12. This idea, later developed by several others (see 
Sadoc and Mosseri |15], Steinhardt et al. jT6], Nelson [T7], Kleman [E]), paved the 
way for a new approach to spatial disorder based on the interplay between order 
and geometry on three-dimensional manifolds of constant Gaussian curvature. 

In-plane order on two-dimensional manifolds has been the subject of much re- 
search since the discovery of the ordered phases and Pp of phospholipid mem- 
branes and is now a rich and mature chapter of condensed matter physics |19] . 
After the seminal work of Nelson and Peliti [20] on shape fluctuations in mem- 
branes with crystalline and hexatic order, much work has been done elucidating 
the intimate relation between in-plane order and the geometry of the underlying 
substrate with many striking results and even more open questions. The funda- 
mental role of topological defects in two-dimensional systems, first elucidated in a 
series of pioneering papers by Berezinskii [21] , Kosterlitz and Thouless \Z'2\ [231 [M] , 
is enhanced in the presence of Gaussian curvature in the underlying medium, and 
leads to structures in the ground state that would normally be highly suppressed 
in flat systems. The goal of this article is to review the most recent developments 
in the study of the ground state properties of two-dimensional order on curved me- 
dia; that is the structure and the mechanics of ordered phases on two-dimensional 
substrates equipped with non-zero Gaussian curvature, in a regime where thermal 
fluctuations are negligible in comparison with other energy scales in the system. 
While focusing on ground states we note that finite temperature physics on curved 
spaces is a rich source of open problems. 

This review is organized as follows. In Sec.[2]we discuss the concept of geometrical 
frustration and review some fundamentals of the elasticity of topological defects 
in flat and curved spaces. The study of crystalline and orientational order on the 
sphere is the basis of most of our current knowledge of order in curved space 
and will be reviewed in Sec. [3j Sec. |4]is dedicated to crystalline order on surfaces 
with variable Gaussian curvature and boundary. The existence of defective ground 
states in toroidal monolayers, with intrinsic crystalline or hexatic order, is a recent 
development in the study of order on curved surfaces and is reviewed in Sec. |5] 
Finally, in Sec. |6] we discuss some current and potential applications of defective 
structures to materials science and nano-engineering. 

1.2. Ordered structures in two-dimensional matter 

Before discussing more technical aspects of the physics of ordered structures on 
curved surfaces we want to recall some salient features of physical systems with 
in-pane order and spatial curvature. 

1.3. Amphiphilic Membranes 

Amphiphilic membranes are thin sheets (50 — 100 A) of amphiphilic molecules im- 
mersed in a fluid and organized in the form of a bilayer (see Fig. [2]). The most 
common constituents of biological membranes are phospholipids consisting of a 
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polar head group and a hydrophobic tail made up of two fatty acyl chains (see Fig. 
[T]). Tails have typical length 14 — 20 carbon atoms and regulate the thickness and 
the stability of the bilayer. The polar head group contains one or more phosphate 
groups -POOH-O-R. Most phospholipid head groups belong to the phosphoglyc- 
erides, which contain glycerol joining the head and the tail. Examples of phospho- 
glycerides include phosphatidylcholine (PC), phosphatidylethanolamine (PE) and 
phosphatidylserine (PS). These are distinguished by the residue R carried by the 
phosphate group (choline: R=-CH2-CH2-N^-(CH3)3 and ethanolamine: R=-CH2- 
CH2-NH2). The fatty acyl chain in biomembranes usually contains an even number 
of carbon atoms. They may be saturated (neighboring C atoms are all connected by 
single bonds) or unsaturated (some neighboring C atoms are connected by double 
bonds). 
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Figure 1. The structure of phosphatidylcholine. 



At low temperature pure phospholipids crystallize and form a bilayer with all the 
tails in the trans configuration and the heads parallel to the bilayer surface and 
firmly linked into a lattice (mostly by hydrogen bonds). At higher temperature 
the order is disrupted and the bilayer becomes fluid. In solution lipid bilayers 
can be found in a number of phases with differing degrees of order among the 
hydrocarbon chains. For phospholipids in the PC family these phases are usually 
called Lq,, Lp and Pp. In the La phase, the tails are liquid and disordered. Lp is 
a gel-like phase in which hydrocarbon tails are ordered and diffusion is severely 
reduced. The Lp La transition is usually referred to as main transition. The 
main transition temperature increases with the length of the molecules due 
to a stronger Van der Waals attraction between adjacent molecules. The degree of 
unsaturation of the hydrocarbon tails also affects Tm- An unsaturated double bond 
can produce a kink in the alkane chain, disrupting the regular periodic structure 
and thus lowering the transition temperature. For bilayers in the PC family, 
ranges form —60 to 80 °C, depending on the tail length and the number of double 
bonds. The order of the hydrocarbon chains also implies a larger thickness of the 
bilayer. The two lamellar phases can be separated by an intermediate "rippled" 
phase Pp in which the bilayer exhibits an undulated structure and almost solid- 
like diffusion properties. Hydrocarbon chains can also appear tilted with respect 
to the bilayer plane. Tilted phases are generally denoted as L'^ and P'p. There are 
in fact several Lp phases characterized by differing amounts of tilt and in-plane 
orientational order [25]. 

Upon changing their concentration, amphiphiles in solution aggregate in a wide 
variety of structures in addition to bilayers. Above the critical micelle concentration 
(which is of the order of 10"'^ mol/£) spherical micelles appear (see Fig. |2]). Their 
formation occurs more readily for single-chain amphiphiles (e.g. monoglycerids) and 
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(C) 

Figure 2. Example of structures formed from self-assembly of amphiphilic molecules, (a) a vesicle, (b) a 
micelle and (c) a bilayer. [Courtesy of Mariana Ruiz Villarreal] . 



is favored by the presence of large head groups. At higher amphiphile concentrations 
spherical micelles are replaced by non-spherical ones and eventually by cylindrical 
rods. Spherical and cylindrical micelles can themselves organize in higher order 
structures such as cubic lattices or hexagonally packed rod piles. More exotic phases 
can be obtained by adding oil to the solution. Once the oil is dispersed in water, 
amphiphiles can form a monolayer across the water-oil interface and self-assemble in 
complex tubular structures known under the common name of plumber's nightmare 
[26]. 

1.3.1. Colloido somes 

The name coUoidosome was coined by Dinsmore et al. to indicate microcapsules 
consisting of a shell of coagulated or partially fused colloidal particles surrounding 
a liquid core [57]. Because of their controllable size, elasticity and permeability, 
colloidosomes have been recognized to form a promising class of "soft devices" 
for the encapsulation and delivery of active ingredients with a variety of potential 
applications for the development of novel drug and vaccine delivery vehicles and 
for the slow release of cosmetic and food supplements. Their major advantage relies 
on the fact that the permeability of the shell depends mainly on the size of the 
gaps between neighboring colloidal particles which can be tuned by controlling 
their size, interactions and degree of fusion. Velev et al. |28| [29] [30] were the first 
to report a method for the preparation of colloidosomes by templating octanol- 
in-water emulsions stabilized by latex particles and subsequently removing the 
octanol core by dissolution in ethanol. Structures of similar architecture have been 
obtained by templating water-in-oil emulsions j3H I32j. Multilayer shells consisting 
of alternating positive and negative polyelectrolytes and / or nanoparticles have also 
been prepared by using layer-by-layer assembly techniques, with the final hollow 
shells being obtained by removal of the central, sacrificial colloidal particles [331 El] • 
Loxley and Vincent [35] developed a new way of preparing polymeric capsules 
with liquid cores based on a phase separation of the polymer within the templated 
emulsion. The colloidosomes produced by Dinsmore et al. [27] were obtained by the 
assembly of latex particles into shells around water-in-oil emulsion drops, followed 
by thermal fusion of the particles in the shell and centrifugal transfer into water 
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Figure 3. Scanning electron microscope image of a coUoidosome from Dinsmore et al. I27| . The coUoi- 
dosome is composed of 0.9 fim diameter polystyrene spheres sintered at 105°C. The close-up on the left 
shows the effect of sintering at the contact points of neighboring spheres. 

through a planar oil- water interface (see Fig. [s]) . 

The coverage of emulsified droplets by colloidal particles takes place by self- 
assembly. The particles dispersed in the fluid spontaneously adsorb on the inter- 
face provided the surface energy between the fluid on the inside and the outside 
of the droplet (cj^o) is larger than the difference of those between the particles 
and the internal fluid (cTp.j) and the particles and the external fluid (o"p_o)- Thus 
Ci,o > |fp,j — Cp^ol- A similar mechanism is used in Pickering emulsions, which are 
stabilized by surface adsorption of colloidal particles. Once adsorbed at the inter- 
face, interacting particles distribute evenly, assuring a full and uniform coverage of 
the droplet. The interactions depend on the type of colloidal particles used as well 
as the liquids. Coated polymethylmethacrylate (PMMA) or polystyrene spheres, 
for instance, acquire a permanent electric dipole moment at the interface between 
the two fluids, probably because of the dissociation of charges on the hydrated 
surface similar to what happens at a water-air interface |36j • The resultant dipolar 
interaction stabilizes the particles and allows full coverage of the droplet. 

The colloidal particles comprising the shell are then locked together to achieve 
the desired permeability and robustness of the colloidosome. Several techniques 
are available to achieve this. By sintering polystyrene particles at a temperature 
slightly above the glass transition {Tg 100 °C), it is possible to achieve a partial 
fusion of neighboring particles at the contact points. This process also allows one 
to control precisely the size of the gaps between particles and therefore the perme- 
ability of the colloidosome. Another method consists of binding the particles with 
a polyelectrolyte of opposite charge which can bridge neighboring particles and 
immobilize them at the interface. Particle locking clearly enhances the toughness 
of the colloidal shell and increases its rupture stress. For sintered polystyrene par- 
ticles the latter can be tuned within the range 1 — 100 MPa. Colloidosomes locked 
with polyelectrolytes are even more deformable and can withstand strains of order 
50 % before rupturing. 

The crystalline arrangement of charged colloids on a hemispherical droplet has 
been recently studied by Irvine and Chaikin who fabricated colloidal suspensions of 
PMMA spheres at the interface between water and cyclohexyl bromide (CHB) [37]. 
Because of the large difference in dielectric constants, ions from the oil strongly 
partition into the water phase [38]. When nearly 100% PMMA particles are dis- 
persed into the oil phase they form a Wigner crystal and a monolayer near the 
interface separated by a zone depleted of particles. These last two features are due 
the ion partitioning that provides the water phase with mobile charges. The net 
charge of water can be controlled through the pH, shifting the hydrolysis and ion 
partitioning equilibria. An electrically neutral water droplet acts as a conductor 
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Figure 4. PMMA colloids sitting at the hemispherical interface between water and cyclohexyl bromide. 
The particles are positively charged and interact via a screened Coulomb interaction with a Debye screening 
length proportional to the concentration of ions in the solvent. [Courtesy of W. Irvine and P. M. Chaikin, 
New York University, New York, NY]. 

and attracts PMMA particles at the interface through an image charge mechanism. 
As a result the particles are barely wet by the water phase and organize in a perfect 
monolayer at (and not across) the interface. 

1.3.2. Viral capsids 

Viral capsids are protein shells that enclose the genetic material of a virus and 
protect it from enzymatic digestion. Capsid proteins are expressed from the DNA 
or RNA genome of the virus and in physiological conditions self-assemble in very 
efficient structures which can withstand high forces (their Young modulus is ~ 2 
MPa) and at the same time effectively disassemble to allow the viral genome to be 
released in the host cell. Most viral capsids have spherical or rod-like shape, but 
less standard shapes, such as conical or toroidal, also occur. 

The crystallographic structures of spherical viruses have been extensively inves- 
tigated and, thanks to the modern techniques of X-ray spectroscopy and cryo- 
transmission electron microscopy, form part of the core knowledge of modern vi- 
rology. In most cases the capsid proteins are grouped in subunits called capsomers, 
oligomers made of either five (pentamer) or six (hexamer) proteins. Spherical 
viruses typically posses icosahedral symmetry with twelve pentamers located at 
the vertices of a regular icosahedron. The number of hexamers that complete the 
capsids is given by 10(T — 1), where T, the triangulation number, takes values 
from a sequence of "magic numbers" (i.e. T = 1, 3, 4, 7. . .) associated with the 
lattice structure of the capsid, as brilliantly explained by Caspar and Klug (CK) 
in a seminal 1962 [39] (the CK construction of icosahedral lattices will be reviewed 
in Sec. [3|. The diameters of spherical viruses span the range from 10 to 100 nm. 
While small capsids are almost perfectly spherical, large viruses, such as the bac- 
teriophage HK97 or the phycodnavirus, typical exhibit a faceted geometry with 
nearly fiat portions separated by ridges and sharp corners corresponding to the 
twelve pentamers. This morphological difference was explained by Lidmar et al 
|40] as a buckling transition resulting from a balance between the stretching energy 
associated with the pentamers in capsomer lattices and the bending elasticity of 
the viral capsid. 

Non-icosahedral capsids of spherocylindrical shape are common among bacterio- 
phages such as some T-even phages as well as the (/)CBK and the 029. In this case 
the capsid has the form of a cylindrical tube composed of a ring of hexamers closed 
at the ends by two half-icosahedral caps. This structure is also found in a vari- 
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Figure 5. (Color online) Two examples of viral capsids. The L-A virus (on the left) and the cowpea mosaic 
virus (CPMV) (on the right). Both have 60 capsomers, but the latter has a marked faceted geometry. Prom 
VIPERdb [46]. 

ant of the T = 7 papovavirus and can be induced in other icosahedral viruses by 
point mutation in the capsid proteins |4H H2] . Of special interest are polymorphic 
viruses, which can appear in either spherical or spherocylindrical conformations. 
Polymorphism has been observed in the polyoma/SV40 animal virus |l3] and the 
cowpea chlorotic mottle virus (CCMV) [S] and, for the latter case, appears to 
be related to the pH and salt concentration of the environment. The human im- 
munodeficiency virus (HIV) also shows broad polymorphism in its capsid shape, 
including cone-like structures in addition to tubular and spherical ones and has 
been well studied in recent years 

1.3.3. Carbon nanotuhes and related materials 

The science of carbon nano-materials has experienced a period of phenomenal 
growth since carbon nanotubes (CNTs) were found by lijima in 1991 [47 and, 
since then, a large number of similar structures, including helix-shaped graphitic 
nanotubes [48,, nanotori [l9], graphitic nanocones [50^ and nanoflowers [51] have 
been reported in the literature. The excitement in this field stems from certain 
exceptional properties that make CNTs potentially useful in many applications 
in nanotechnology, electronics, optics and other fields of materials science. They 
exhibit extraordinary strength, unique electrical properties and are efficient con- 
ductors of heat. Several methods for the preparation of graphitic nanomaterials 
have been developed, including laser pyrolysis, arc discharge, and electron irradia- 
tion. Recently, metal-catalyzed methods have also been used to synthesize carbon 
nanomaterials. 

Most single wall nanotubes (SWNT) have diameter of about 1 nm, while the 
length is often of the order of microns. The lattice structure of a SWNT can be 
obtained from that of a graphene plane by assigning a pair of indices (n, m) which 
specify how the graphene lattice is rolled up into a seamless cylinder (see Sec. 
|5]). n = nanotubes are referred to as "zigzag", while n = m tubules are called 
"armchair" . The generic name "chiral" is used otherwise |52] . In terms of tensile 
deformations, SWNT are the stiffest materials known with a Young modulus in 
the range 1-5 TPa and a tensile strength of 13-53 GPa. This strength results from 
the covalent sp^ bonds formed between the individual carbon atoms. 

The lattice structure of a nanotube strongly affects its electrical properties be- 
cause of the interplay between the unique electronic structure of graphene and the 
tubular geometry. For a given (n, m) structure, a nanotube can be a conductor 
(if n = m), a small-gap semiconductor (if n — m is a multiple of 3) or a standard 
semiconductor (otherwise). Thus all armchair (n = m) nanotubes are metallic, and 
nanotubes with {n,m) equal (5,0), (6,4), (9,1), etc. are semiconducting. 

Topological defects may occur on the side-wall of carbon nanotubes and related 
materials in the form of atomic vacancies, 5 — 7 dislocations and Stone- Wales defects 
|53] (i.e. quadrupoles consisting of two pairs of 5-membered and 7-membered rings) 
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and are believed to significantly change the mechanical and transport properties 
of carbon nanotubes. The latter, in particular, were suggested to serve as possible 
nucleation centers for the formation of dislocations in the original ideal graphene 
network and constitute the onset of possible plastic deformations [53^ . Electronic 
transport is affected by the presence of defects by lowering the conductivity through 
the defective region of the tube. Defect formation in armchair-type tubes (which 
can conduct electricity) can cause the region surrounding that defect to become 
semiconducting. Furthermore, single monoatomic vacancies induce magnetic prop- 
erties. Phonon scattering by defective regions heavily affects the thermal properties 
of carbon nanotubes and leads to an overall reduction of the thermal conductivity. 
Phonon transport simulations indicate that substitutional defects such as nitrogen 
or boron will primarily lead to scattering of high frequency optical phonons. Larger 
scale defects, however, such as Stone Wales defects, cause phonon scattering over 
a wide range of frequencies, leading to a greater reduction in thermal conductivity 
|55] . On the other hand, defective regions appear to be natural places for chemical 
functionalization. Numerical simulations have shown how the presence of Stone- 
Wales defects considerably enhances the adsorption of carboxyl groups (COOH) 
which can then bind to molecules with amide and ester bonds |5B] . 

2. Interacting topological defects in curved media 
2.1. Geometrical frustration 

The notion of geometrical frustration was introduced to describe situations where 
certain types of local order, favoured by physical interactions, cannot propagate 
throughout a system [57] . The expression was used for the first time by Toulouse in 
1977 [58j to describe certain particular magnetic systems with nearest-neighbours 
interactions which cannot all be satisfied simultaneously. A textbook example of 
frustration in magnetic models is represented by a system of Ising spins on a 
triangular lattice with antiferromagnetic bonds: while a perfect antiferromagnetic 
alignment would minimize all terms in the Ising Hamiltonian, such an alignment 
is not allowed by the topology of the underlying lattice so that for any triangular 
plaquette there is always at least one unsatisfied bond (see Fig. |6]) . This concept 
can be extended naturally to any system where interactions impose a local order, 
but the most favoured local configuration is geometrically incompatible with the 
structure of the embedding space. 

Two-dimensional manifolds equipped with some microscopic field for which a 
notion of local order can be defined unambiguously provide a paradigm for systems 
exhibiting geometrical frustration. Consider for instance an assembly of identical 
particles interacting with a spherically symmetric pair potential Vij = V{\xi — Xj\)^ 
with Xi the position vector of the ith particle in a suitable coordinate system. 
In flat two-dimensional space, particles almost always pack in triangular lattices, 
unless the interaction potential is carefully tuned to select some some other lattice 
topology. Endowing the medium with a non-planar topology introduces frustration 
in the sense that the energetically favoured 6— fold orientational order can no longer 
be established everywhere in the system. Such geometrical frustration arises at the 
microscopic level from the celebrated Euler theorem of topology which relates the 
number of vertices V , edges E and faces F of any tessellation of a 2— manifold M: 

V-E + F = x, (1) 
where x is the Euler characteristic of M. If M is an orientable closed surface, one 
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Figure 6. (Color online) Examples of geometrical frustration. (Left) Ising antiferromagnet on a triangular 
lattice. Because of the topology of the underlying lattice, no arrangement of the three spins can minimize 
the Hamiltonian H = —J^li^ij) SiSj. (Center) Vector field of a sphere. As a consequence of the Poincare- 
Hopf theorem, a vector field must vanish at least at two points, corresponding to the north and south 
pole of the sphere in this example. (Right) Triangulation of the sphere. As required by the Euler theorem, 
any triangulation of the sphere must contain some vertices of coordination number different from six (i.e. 
twelve 5— fold vertices in this case). 

can show that x is determined uniquely by an integer g > 0, cahed the genus of 
M, which represents the number of "handles" of M; namely x = 2(1 — g). Two 
orientable closed surfaces with the same genus (thus the same Euler characteristic) 
are homeomorphic: they can be mapped into one another without changing their 
topological properties. In a surface with boundary, the Euler characteristic is given 
hy X = 2(1 — g) — h, where h is the number of boundaries or "holes" of M. Thus 
a sphere, which has no handles nor boundary (g = h = 0) has x = 2, while the 
embedded torus {g = 1 and h = 0) has x = 0. A disk, on the other hand, has x = 1 
{g = and h = 1) and is topologically equivalent to a sphere with one hole. In 
Sees. [3](5] we discuss ordered structures on three important topologies: the sphere, 
the disk and the torus. 

In the case of two-dimensional crystals with 6— fold local order Eq. ([T| can be 
rephrased in a form that is particularly useful in describing the presence of defects 
in the lowest energy state by defining a topological charge as the departure from 
the ideal coordination number of a planar triangular lattice: qi = 6 — Cj, with q the 
coordination number of the ith vertex. Now, consider a tessellation in which each 
face is an n— sided polygon and let k faces meet at each vertex. Since each edge is 
shared between two faces and links two vertices it follows that: 

nF = 2E = Y^ kVk , 
k 

where Vk is the number of vertices of degree k. For a triangulation n = 3. From 
Eq. ([T]) it follows then: 

V 

Q = ^q, = 6x. (2) 

i=l 

In the case of a sphere, with x = 2, Eq. (|2| implies any triangulation contains 
defective sites such that the total topological charge of the lattice is Q = 12. This 
can be achieved, for example, by incorporating twelve 5— fold disclinations (with 
g = 1) in a network of 6— fold coordinated sites like in a common soccer ball. 
These twelve disclinations, which would be suppressed in the lowest energy state 
of a planar crystal, result from the geometrical frustration associated with the 
topology of the sphere. 

Eq. (U is a special case of geometrical frustration on 2-manifolds where local 
orientations are defined modulo tt/3. More generally one can consider a p— atic 
director field for which local orientations are defined modulo 2tt/p. The topological 
charge of a disclination in this case is q = AO / ^ , where is the angle the director 
rotates in one counter-clockwise circuit of any closed contour enclosing the defect. 
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Figure 7. (Color online) Examples of 5— fold (left) and 7— fold (right) disclinations on a triangular lattice 
in regions of positive and negative Gaussian curvature respectively. 



where is the total number of defects. The ratio k = q/p is commonly referred to 
as the winding number of the disclination. In the case of a simple vector field, for 
instance, p = 1 and Eq. Q corresponds to the well known Poincare-Hopf theorem 
according to which the sum of the indices of all the isolated zeros of a vector field 
on a oriented differentiable manifold M is equal to the Euler characteristic x of 
M. Thus for a sphere a vector field must have at least one sink and one source, 
each having topological charge one, while on a torus (x = 0) a vector field can 
be defect free. A nematic director n, on the other hand, has p = 2 (i.e. physical 
configurations are invariant under inversions n — > — n). Thus disclinations with ±1 
topological charge correspond to configurations where the director rotates itTr in 
one circuit enclosing the defect. These elementary disclinations have semi-integer 
winding number k = ±1/2. As a consequence of Eq. Q the total topological 
charge of a nematic texture on the sphere is Q = 4, corresponding for instance to 
the typical baseball texture consisting of four q = 1 {k = 1/2) disclinations located 
at the vertices of a regular tetrahedron [591 EO] • 

2.2. Mathematical preliminaries and notation 

The equilibrium structure of two-dimensional locally ordered systems on curved 
substrates depends crucially on the existence and arrangement of the defects. The 
seminal work of Kosterlitz, Thouless, Halperin, Nelson and Young (KTHNY) on 
defect-mediated melting in two dimensions makes it clear that it is useful to em- 
ploy a theoretical framework where the fundamental objects in the system are the 
defects themselves, while treating the microscopic constituents within a continuum 
elastic theory. This approach has the advantage of far fewer degrees of freedom 
than a direct treatment of the microscopic interactions and allows one to explore 
the origin of the emergent symmetry observed in non-Euclidean ordered structures 
as a result of the interplay between defects and geometry. This interplay is one 
of the fundamental hallmarks of order on two-dimensional manifolds and leads 
to universal features observed in systems as different as viral capsids and carbon 
macromolecules. In this section we will briefly review some concepts of differential 
geometry, mostly to establish notation. In the next section we will review some fun- 
damentals of the elasticity of defects in two dimensions. The coupling mechanism 
between curvature and defects will be introduced in Sec. 12.41 
Points on a two-dimensional surface S embedded in are specified by a 



Eq. ([2]) becomes then: 



N 




(3) 



1=1 
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three-dimensional vector R{x) as a function of a two-dimensional coordinate 
X = (x^, a;^). For each point of S we define three vectors: 

Qi = diR 1 = 1,2 (4) 

and n = ^ — , (5) 
Iffi X ff2| 

where di = d/dx"^. The vectors Qi belong to Tr5, the tangent space of 5 at 
whereas n is a normal vector. Note that while n is a unit vector, gi are generally 
not of unit length. The metric, or first fundamental form, of S is defined as: 

ds^ = gijdx^dx^ (6) 

where gij is the metric tensor: 

Qij = Qi ■ Qj (7) 
The dual tensor is denoted as g'^^ and is such that: 

gikg^^ = Sf , 

with 6j the Kroncckcr symbol. This allows us to introduce contravariant tangent- 
plane vectors gr* = g^-^gj, satisfying gi - g^ ~ H - -^^y vector v on the tangent plane 
can be expressed as a linear combination of basis vectors gi oi g^: v = v^gi = v^g"^, 
where Vi = gijvK The extrinsic curvature of the surface S is encoded in the second 
fundamental form bij (also known as the extrinsic curvature tensor): 

bij = -Qi ■ djU = n ■ djgi . (8) 

The eigendirections of bij at a given point correspond to the principal curvature di- 
rections of S at that point and the associated eigenvalues ki and H2 are the extremal 
(or principal) curvatures. The mean curvature H and the Gaussian curvature K 
are defined as the sum and the product of the principal curvatures: 

2H = Ki + K2 = gijb'^ (9) 
K = K1K2 = l^^eHijbki (10) 

where e^^ is the dual of the Levi-Civita tensor, whose components are given by: 

en = €22 = 
ei2 = -£21 = \/g 

where g = det^jj. The contravariant form is e*-' = eij/g and satisfies Sike'^ = Sj. 
Since CijV^v^ = 0, where is any contravariant vector, it follows that the vector 
eijV^ is perpendicular to . Thus inner multiplication by e^j rotates a vector by 
7r/2. The covariant derivative of a vector field v in the zth coordinate direction is 
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defined as usual by: 

Viv'' = div'' + T'^jV^ (11a) 
^iVk = diVk-Tljjjj (lib) 
where V^^ is the Christoffel symbol: 

rf,- = ^5^' {dm + di9lj - di9ij) . (12) 

Both the metric and the Levi-Civita tensor are invariant under parallel transport. 
This translates into: 

^k9ij = ^k9'' = VkCij = Vke'^ = . 

Much of the elastic theory of defects, either in fiat or curved systems, relies on the 
calculation of the Green function of the Laplace operator. On a generic 2— manifold 
the latter obeys: 

/^GL{x,y) = 5{x,y), (13) 
where A is the Laplace-Beltrami operator: 

A = ^d,{^g^^d,) (14) 

and 5 the delta-function: 

(15) 

V 9 

The Stokes theorem is frequently invoked when calculating elastic energies of de- 
fects. In covariant form it can be stated as follows: given a vector field w on a 
2— manifold M with boundary dM, the following identity holds: 

/ dx^Vk= [ (fxe'^ViVj. (16) 

JdM Jm 

For consistency we will adopt covariant notation throughout this review. When 
discussing planar systems, in particular, we have: bij = H = K = and the 
elements of the metric tensor are given in Cartesian coordinates by gxx = 9yy = 1 ; 
9xy = 9yx = and in polar coordinates by i?rr = 1, 9h = r"^-, 9r<l> = 9(l>r = 0. 

The definition of oricntational order on a surface clearly requires a non- 
ambiguous notion of angular distance between vectors on the same tangent plane. 
This is traditionally achieved by introducing a pair of orthonormal vectors Bq, 
{a = 1, 2) called vielbein (note that the canonical coordinate vectors gi are gener- 
ally neither orthonormal nor orthogonal) so that: 

ea • 6/3 = (ea)i(e^)* = Sa/3 , (ea)i(ea)i = 9ij , (17) 



where Sap is the usual Kronecker symbol in the indices a and (3. A vector field 
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V = v^Qi can be expressed alternatively in the basis e^: 

V = v"ea v"' = v''{ea)i ■ 

Since the coordinates v"' are locally Cartesian, dap = S'^^ and there is no distinction 
between upper and lower Greek indices: Va = v"'. Vielbein are constructed to be 
invariant under parallel transport, and thus Vi{ea)j = and: 

ViVa = {ea)jViV^ = diVa + ^iafiVp , 

where ^iaj3 is the so called spin connection, and is given by: 



iaf3 — ' 

Taking the derivative of the left equation in ( |17| ) one immediately sees: 

{ea)''di{ef3)k = -{ef3fdi{ea)k , 

from which it follows that Oj^^ = — rjj^^- In two dimensions this permits the 
parametrization of the spin-connection by 8' single covariant vector Q, (i.e. 
the antisymmetry under exchange of the Greek indices leaves only — l)/2 = 2 
independent components in d = 2 dimensions). Thus 



(18) 



with 



— (5o(^^ the antisymmetric symbol. 



to/3 — ^a'-'p "fj'^a 

It is well known how the curvature of a manifold manifests itself when a vector 
is parallel transported around a closed path. Taking an infinitesimal square loop of 
sides dx and dy, the parallel transported vector v' differs from the original vector 
V by an amount: 



RijvUx'dy^ , 



where R^- is the Riemann tensor given in two-dimensions by: 



ijk 



KeU 



jk 



(19) 



(20) 



In the orthonormal basis e^, Eq. (19) reads: 



v'a-Va = Rijapvpdx'dy^ , 
where Rija/s is the curvature tensor associated with the spin-connection rjja/j: 



In two dimensions, using Eq. (20), one can write: 

Rijai3 = {ea)^{eii)^ Rijki = {ea)^{epfeijekiK = eais^ijK , 



which, combined with Eqs. (21) and (18), yields 



Rijal3 — ^cxpidi^j — dj^i 



(22) 
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Figure 8. Examples of topological defects in a vector field, (a), (b) and (c) have topological charge q - 
(d) q = —I and (e) and (f) have charge q = 3 and q = —3 respectively. 



which also imphes: 

V X n = e'^Vi^j = K . (23) 



2.3. Elasticity of defects in the plane 

The Xy— model is the simplest setting where particle-like objects emerge from 
a purely continuum theory [22l |23l |2l] . The order parameter is the angular field 
6 S [0,27r], which may represent the orientation of two-dimensional vectors S = 
S{cosO,sm.O) or the phase of a complex field = |V'|e*^- The interaction that 
tends to minimize the spatial variations of the order parameter results from the 
continuum free energy: 

Fei = \KA j cPx\Ve{x)\^ . (24) 

Despite its simplicity this model successfully describes several aspects of the physics 
of vortices in superfiuid '^He or thin superconducting films, where the angle 6 is 



identified with the phase of the collective wave function. Eq. (24) also describes 
nematic liquid crystals in the limit of equal splay and bending moduli. This ap- 
proximation, however, favors q = ±2 {k = ±1) disclinations rather than the more 
natural q = ±1 (A; = ±1/2) disclinations of nematics and is not particularly suit- 
able to describe the ground state. At T > 0, thermal fluctuations drive the two 
elastic constants to the same value at long wavelengths, so that there is a unique 
Kosterlitz-Thouless transition temperature [61]. As discussed by Deem [62], the 



essential effect of unequal elastic constants is to create a distinct long-range contri- 



bution to the core energy of each defect. More generally Eq. (24) can be considered 



as the simplest phenomenological free energy describing p— atic order (assuming 9 



defined modulo 27r/p). Calling v = V0, the minimization of the free energy (24) 
leads to the equilibrium condition: 



= Viv' = . 



(25) 
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In the presence of a number of disclinations of topological charge qa, 9 changes 
by (27r/p) (^Q, in one circuit around any contour enclosing a total topological 
charge la- 

dx' Vi = — ^Qa- (26) 



Using Stokes' theorem, Eq. (26) can be translated into the requirement: 

e'iV.Vj = e'^ViVjO = r){x) , (27) 

where: 

ri{x) = — y^qad{x,Xa) 

a 

is the topological charge density. Using standard manipulations (see for example 
[63]) a vector field v satisfying Eq. (25) with the constraint (26) can be found in 
the form: 

Vi{x) = -4V, I cfy Gl{x, y)ri{y) , (28) 

where GL{x,y) is the Laplacian Green function. Using Eq. (28) in Eq. (24) leads 
to the well known expression: 

Fei = -^Ka J <fxcfiyGL{x,y)7]{x)ij{y). (29) 

Like charged particles, disclinations interact via a Coulomb potential, which in two 
dimensions is proportional to the logarithm of the distance in the plane. As noted, 
this framework is valid for both nematic liquid crystals in the one elastic constant 
approximation and superfluids. An important difference between these two systems 
lies in the choice of the boundary condition for the field 6. Nematogens are typically 
forced to be normal to boundary of the substrate and this implies a constraint for 
6, while there is no such constraint for ^He films since the wave function is defined 
in a different space from that to which the superfluid is confined. This difference 
is crucial on a curved substrate (see Ref. [M] for a detailed review of this topic). 

Eq. (29) represents the elastic energy associated with the distortion introduced 
by defects in the far field, where the elastic variables change slowly in space. This 
expression breaks down in the neighborhood (or core) of a defect, where the order 
parameter is destroyed and the actual energetic contribution depends on micro- 



scopic details. In order to describe a defective system at any length scale, Eq. (29) 
must be corrected by adding a core energy Fc representing the energetic contribu- 
tion within the core of a defect, where standard elasticity breaks down. A detailed 
calculation of the core energy requires some microscopic model and is usually quite 
complicated. Nonetheless its order of magnitude can be estimated by writing: 

Fc = vraVc , (30) 

where a is the core radius, corresponding to the short distance cut-off of the elastic 
theory, and fc is some unknown energy density independent of a. fc can then be 
estimated by minimizing the total energy of the system with respect to a. For a 
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planar system with a single disclination of winding numer /c, the total energy can 



be easily calculated from Eqs. (29) and (30): 



(31) 



Minimizing Eq. (31) with respect to a, one finds fc = i^Afc^/(2a^) from which: 

1 



(32) 



The quantity ec is the energy needed to increase the number of defects by one, 
regardless of its location. 

The elasticity of dislocations and disclinations in solids resembles in many re- 
spects that of vortex lines in the Xl"— model. In two dimensional elasticity a pure 
in-plane deformation is encoded in a displacement field u*, i = 1, 2, which maps 
any point in the system x to: 



X = X + U Qi 



(33) 



where gi is a suitable basis in the coordinates of the undeformed system. If there 
are no defects, the displacement field is a single- valued mapping of the plane into 
itself. Topological defects introduce an incompatibility in the displacement field, 
in the sense that is not a single-valued mapping anymore. The elastic stress in 
the region surrounding a defect appears at the macroscopic level from the Hooke's 
law of elasticity: 



l + v 
Y 



Y 



9ij<^k 



(34) 



where Y and v are the two-dimensional Young modulus and Poisson ratio re- 
spectively and Gij is the stress tensor. Eq. (34) can be obtained, for example, by 
minimizing the elastic energy: 



(35) 



where A and /i are the Lame coefficients in two dimensions: 

Yv Y 



A 



2{l + u) 



In the absence of body forces, the force balance equation requires aij to be diver- 
gence free: 



Vja'^ = . 



(36) 



The strain tensor mj represents the variation in the first fundamental form of the 
surface due to the deformation field (33), namely: 



2ui 



gij{x + u) - gij{x) = ViUj + VjUi + 0{u^) . 



(37) 
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In the presence of a dislocation line L the function u becomes multivalued so that, 
while traversing any closed counterclockwise loop C containing L: 



du 



(38) 



c 



where h is the Burgers vector representing the amount by which the image of a 



closed loop under the mapping (33) fails to close in presence of a dislocation line 



|65|I66]. If L is an isolated straight dislocation with origin at the point a;o, Eq. (38) 
implies: 



(39) 



Generally speaking, maps such as that of Eq. (33), induce variations in both the 



length and the orientation of an infinitesimal distance vector dx in the deformed 
medium. This statement can be clarified by writing: 



dx[ — dxi = (VjUi) dx^ 



{uij -LOij) dx^ 



where: 



^ij = 2 (VjUj - VjUi) 



(40) 



is the infinitesimal rotation tensor induced by the deformation Eq. (33). In two 
dimensions the latter can be conveniently written in the form: 



UJi 



9 



(41) 



with: 



G 



(42) 



Since the application of dj to a vector rotates the vector by 7r/2 clockwise, the 



coefficient in Eq. ( 42 ) can be interpreted as the average rotation of an infinitesi- 



mal line element under a distortion of the form Eq. ( 33 ) . In presence of an isolated 



disclination the field & becomes multi-valued. If s is the angular deficit associated 
with the disclination, integrating along a circuit enclosing the disclination core 
yields: 



de = s , 



(43) 



c 



which can be rephrased as a statement about commutativity of partial derivatives: 



e'^ViVjO = s6{x,xo) 



(44) 



where xq is the position of the disclination core. In a lattice with n— fold rotational 
symmetry, s = {2TT/n)q, where q is the topological charge introduced above. A 



comparison between Eq. I\27h and Eq. (44) clarifies the common identification of Q 



with the bond angle field of a crystal [67^ . 



Eqs. (39) and (44) provide a set of relations between the fundamental features of 



topological defects in crystals (i.e. the Burgers vector b and the topological charge 
q) and the partial derivatives of the displacement field u. These relations can be 
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used to derive an equation for the elastic stress arising in a system as a consequence 
of the defects. A simple and elegant way to achieve this task in the context of planar 
elasticity relies on the parametrization of the stress tensor via a single scalar field 
X known as the Airy stress function (see for example p)8!). Taking advantage of 
the commutativity of partial derivatives in Euclidean space, one can write 



a ■ 



(45) 



so that the force balance condition (36) is automatically satisfied. Applying the 
operator e*'^e-''VfcV/ to both sides of Eq. (34) and using Eq. (45) one finds the 

(46) 



following fourth order Poisson-like problem: 

/\''x{x) = Yr]{x) , 

where is the biharmonic operator and rj is the defect charge density: 

2tt 



n 



where Sq, denotes the topological charge of a disclination at x^ and 6^ the Burgers 
vector of a dislocation ai x^. 

Let's now turn our attention to the case of a planar crystal with local 6— fold 
orientational symmetry populated with disclinations of density: 



r]{x) 



N 

-y 



qaS{x,Xa) . 



(47) 



Assuming a free boundary (i.e. all the components of the stress tensor are zero 
along the boundary and thus x = = 0); the stress function x can be expressed 
in the Green form: 



Y 



(fyG2L{x,y)ir]{y) , 



(48) 



where G2l{x, y) is the biharmonic Green function. The elastic energy of the system 
is given by Eq. (pSj) and (48) in the form: 



^ei = ^Y j (fx(fyG2L{x,y)r]{x)r]{y). 



(49) 



Eqns. (24) and (49) are the fundamental equations in the elastic theory of defects 
in two-dimensional planar p— atics and crystals. In the next section we will see how 
a non-zero Gaussian curvature in the underlying medium affects these energies by 
effectively screening the topological charge of the defects. 



2.4. Coupling mechanisms between curvature and defects 

The elasticity of topological defects on curved substrates equipped with local ori- 
entational order was first considered by Nelson and Peliti in the context of hex- 
atic membranes [20]. It is now standard knowledge in the statistical mechanics 
of membranes that long range forces appearing as a consequence of local orienta- 
tional order crucially affect the behavior of membranes at finite temperature. The 
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stiffness associated witli orientational correlations leads to an enhancement of the 
bending rigidity that counteracts the thermal softening that occurs in fluid mem- 
branes. Furthermore, for planar membranes, the stabilizing effect induced by the 
orientational stiffness opposes the entropy-driven tendency to crumple^ causing a 
transition between a flat and crumpled phase at T > 0. In this review we focus on 
the ground state properties of ordered structures on curved surfaces and we refer 
the reader to the literature for a discussion of finite temperature physics [19j. 

In this section we will show how an underlying non-zero Gaussian curvature cou- 
ples with the defects by screening their topological charge. As a result, topological 
defects, whose existence would be energetically suppressed if the same system was 
lying on a flat substrate, instead proliferate. The universality of such a curvature 
screening mechanism is remarkable in the sense that it occurs in a conceptually 
identical fashion in both p— atics and solids, although the kernel of the elastic inter- 
actions between defects differs. As a consequence of curvature screening topological 
defects organize themselves on a rigid surface so as to match the Gaussian curva- 
ture of the substrate. On the other hand, if the geometry of the substrate is allowed 
to change (for instance by lowering the bending rigidity) the system eventually en- 
ters a regime where the substrate itself changes shape in order to accommodate 
some preferential in-plane order. The latter is the fundamental mechanism behind 
the buckling of crystalline membranes [HS] and is believed to be the origin of the 
polyhedral geometry of large spherical viruses such as bacteriophages [7D] . 

In the case of manifolds equipped with a pure rotational degree of freedom like 
p— atics, the curvature affects the elastic energy through the connection which 
determine how vectors change when parallel transported. This statement can be 



clarified by rewriting the elastic energy (24) in the form: 



Fei = ^Ka j (fx Viin^V'mj , (50) 

where m is a p— atic director field which can be conveniently expressed in a local 
orthonormal frame (a = 1 2): 

m = cos 9 Bi + sin 9 62 , 

with 6 defined modulo 27r/p. Since dinia = —ectprapdi9^ using the properties of 
vielbein outlined in Sec. 12.21 we have: 



ViTTij = {ea)j{dima + Vtiapmp) = -eapm,3{ea)jidie - (51) 



The elastic energy (50) thus becomes 



Fei = ^Ka j d'xg'\di9 - ni){dj9 - n,) . (52) 
Using Eq. (23), the vector field can be expressed as: 

= -4^J J Gl{x, y)K{y) , (53) 
where GL{x,y) is again the Green function of the Laplace-Beltrami operat or (|14[ ). 



Substituting Eqns. (53) and (28) in the expression for the elastic energy (52) we 
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obtain: 



Fei = -\ka j d^xd^yGL{x,yMx) - K{x)My) - K{y)] . (54) 
As anticipated, the Gaussian curvature of the underlying substrate couples with the 



defects by screening their topological charge. Eq (54) is identical to the Coulomb 
energy of a multi-component plasma of charge density ry in a background of charge 
density —K. Like particles in a plasma, we can expect defects to screen each other's 
charge, thus forming clusters of zero net charge, and matching the charge distribu- 
tion of the surrounding background. This implies that disclinations will be attracted 
to regions of like-sign Gaussian curvature. 

Crystalline surfaces (i.e. non-Euclidean crystals) differ from manifolds equipped 
with a p— atic director field because bonds, whose local orientation is encoded 
in the field 9 in p— atics, can now be compressed and sheared to relieve part of 
the elastic stress due to the geometrical frustration provided by the embedding 
manifold. A first attempt at describing the elasticity of defects on a two-dimensional 
curved crystal was made by Dodgson ^71j, who studied the ground state of the 
Abrikosov fulx lattice in a model thin-film superconductor on a sphere (subject to 
a field radiating from a magnetic monopole at the center) and found evidence for 
twelve 5— fold disclinations at the vertices of an icosahedron in a otherwise 6— fold 
coordinated environment (see Sec. [s] for a detailed review of spherical crystals). 
Later, Dodgson and Moore [72] proposed adding dislocations to the ground state of 
a sufficiently large spherical vortex crystal to screen out the strain introduced by the 
twelve, topologically required, disclinations. A general framework to describe the 
elasticity of defects in non-Euclidean crystals was proposed by Bowick, Nelson and 
Travesset (BNT) in 2000 [73]. The BNT model, obtained from the covariantization 



of the elastic energy ( 49 ) , relies on the following expression for the elastic energy of 
a collection of disclinations in a triangular lattice of underlying Gaussian curvature 
K: 

F,i = ^Y I d^xd^yG2L{x,yMx) - K{x)My) - K{y)] , (55) 



where "q^x) is the topological charge density (47) and G2l{x, y) the Green function 
of the covariant biharmonic operator on the manifold. The origin of the coupling 
between topological charge and curvature, in this case, is rooted in a profound 
result of discrete geometry originally due to Descartes which can be considered the 
oldest ancestor of the Gauss-Bonnet theorem of differential geometry. Let P be a 
convex polyhedron and define the angular deficit of a vertex v of P as: 



c{v) 

k{v) = 27r-^ai{v), (56) 

i=l 

where ai{v) with i £ [l,c{v)] are the angles formed by all the faces meeting at v. 
Descartes' theorem states that the sum of the angular deficits of a convex polyhe- 
dron is equal to Air. More generally: 

^ k{v) = 27rx , (57) 

which is exactly the Gauss-Bonnet theorem for the case in which the Gaussian cur- 
vature is concentrated in a finite number of points (vertices) rather than smoothly 
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Figure 9. (Color online) Parallel transport of a vector around a corner of a cube. As a consequence of the 
discrete Gaussian curvature, the vector is rotated by 7r/2 after parallel transport. 

distributed across the whole surface. The deficit angle k{v) is thus the discrete 
analog of the Gaussian curvature. This analogy is not limited exclusively to Eq. 
(56). Consider, for example, the corner of a cube and imagine parallel transporting 
a vector v around a closed loop encircling the corner (see Fig. [9|. After parallel 
transport, the vector has rotated by k{v) = 7r/2 with respect to its original orienta- 
tion. Analogously, a vector parallel transported around a closed loop on a surface 
rotates by an angle 



where the integral is over the portion of the surface enclosed by the loop. Now, on 
a triangulated surface a{v) ~ vr/S and: 



TT 



k{v) = 2tt — — c{v) = — [6 — c{v)] 



Thus the source term rj — K figuring in Eq. (55) corresponds to the difference 
between the pre-existing Gaussian curvature of the manifold with local 6— fold 
orientational order and the additional Gaussian curvature induced by a disclination 
of topological charge q. A more formal discussion can be found in Ref. [71] . As in the 
case of a purely rotational degree of freedom, disclinations in non-Euclidean crystals 
organize so as to match the Gaussian curvature of the underlying medium while 
maximizing their reciprocal distance as a consequence of the repulsive interaction 
between like-sign defects. In the next three sections we will see how the formalism 
outlined here applies to three physically relevant examples of manifolds: the sphere, 
the topological disk and the torus. 



3. Crystalline and nematic order on the sphere 
3.1. Introduction 

Experimental realizations of spherical crystals are found in emulsions of two im- 
miscible fluids, such as oil and water, stabilized against droplet coalescence by 
coating droplets of one phase (say water) with small colloidal particles [75l [76] . 
The colloidal "armor plating" of these Pickering emulsions also plays a role in col- 
loidosomes, colloid-coated lipid bilayer vesicles used for encapsulation and delivery 
of flavors, fragrances and drugs [27j. Identical micron-sized particles tend to crys- 
tallize under typical experimental conditions. The defect structure of crystalline 
ground states on a sphere is important in this context since defects influence the 
strength of the colloidal armor. 
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Defects play an essential role in describing crystalline particle packings on the 
sphere. At least twelve particles with 5-fold coordination (i.e., 12 disclinations) 
are required for topological reasons, and like the 5-fold rings in carbon fullerenes, 
one might expect that the energy would be minimized if the disclination positions 
approximated the vertices of a regular icosahedron. This expectation, which also 
plays a role in geodesic domes and in the protein capsomere configurations of 
spherical virus shells [39| [77] , is nevertheless violated when the shells are sufficiently 
large and disclination buckling p[Q] out of the spherical environment is suppressed 
by surface tension. 

To understand the structure of spherical crystals, we must find the ground state of 
M particles distributed on a surface S and interacting with, say, pairwise repulsive 
potentials. For simplicity consider power law potentials V{r) = e^/r'^, where e 
is a generalized electric charge. The energy of M particles at positions r{i), with 
(■ = {h^h) e becomes 2Eq = Y^t^e /\r{tj-r{i')\^ . This covers most of the 

interesting physical systems. The case 7 = 1 corresponds to the pure 3D Coulomb 
potential, appropriate for charged colloids. This is known as the Thomson problem 
in the physics and mathematics literature and has been widely studied for over 100 
years [751 ES HO]- The case 7 = 8 corresponds to a dipole interaction, appropriate 
for neutral colloids at the interface between two liquids. The different dielectric 
constants of the two liquids leads to an asymmetric distribution of charge on the 
colloids and a net dipole moment. The case 7 = 12 is the repulsive part of the 
Lennard-Jones potential and is the important piece of the interaction for driving 
crystallization. Finally the limit 7 —> 00 is an extreme short range interaction that 
selects the pair of particles with the minimum distance. The equivalent packing of 
spherical caps was introduced in 1930 by the biologist Tammes [8T] and connects 
the present problem to the rich physics and mathematics of packing [82] |83l . 

The packing of particles on surfaces of nontrivial topology is relevant to a broad 
range of physical, chemical and biological systems in addition to its intrinsic math- 
ematical interest. An almost literal realization of the Thomson problem is provided 
by multi-electron bubbles [85] . Electrons trapped on the surface of liquid helium 
have long been used to investigate two dimensional melting j86[ 157] . Multi-electron 
bubbles result when a large number of electrons (10^—10^) at the helium interface 
subduct in response to an increase in the anode potential and coat the inside wall of 
a helium vapor sphere of radius 10— lOO^m. Typical electron spacings, both at the 
interface and on the sphere, are of order 2000 A, so the physics is entirely classical, 
in contrast to the quantum problem of electron shells which originally motivated 
J.J. Thomson [75]. Information about electron configurations on these bubbles can, 
in principle, be inferred from studying capillary wave excitations |88| 189] . Similar 
electron configurations should arise on the surface of liquid metal drops confined 
in Paul traps [90'. 

A Thomson-like problem also arises in determining the arrangements of the pro- 
tein subunits which comprise the shells of spherical viruses [391 HOI HSj [771 [9l] . Here, 
the "particles" are clusters of protein subunits arranged on a shell. This problem of 
protein arrangements was solved in a beautiful paper by Caspar and Klug |i39j for 
intermediate values of R/a, where R is the sphere radius and a is the mean particle 
spacing. Caspar and Klug constructed icosadeltahedral particle packings charac- 
terized by integers m and n, which provide regular tessellations of magic numbers 
M = 10(n^ -|- nm -\- w?) + 2 particles on the sphere. The Caspar-Klug tessellations 
provide an excellent starting point for finding low energy particle configurations 
on the sphere for intermediate values of Af ~ 87r(ii/a)^/\/3- Particle numbers M 
that are not magic numbers can be accommodated by introducing vacancies or 
interstitials into these regular packings. 
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Other realizations of Thomson-like problems include regular arrangements of 
colloidal particles in colloidosome cages [92l |93l [Ml ES] proposed for protection of 
cells or drug-containing vesicles [27] and fuUerene patterns of spherical carbon 
atoms [96 . An example with long range (logarithmic) interactions is provided 
by the Abrikosov lattice of vortices which would form at low temperatures in a 
superconducting metal shell with a large monopole (solenoid tip) at the center 

The problem of the best possible packing on spheres also has applications in the 
micro-patterning of spherical particles [97], relevant for photonic crystals, and in 
understanding the structure of Clathrin cages, responsible for the vesicular trans- 
port of cargo in cells [98] . Crystalline domains covering a fraction of the sphere are 
also of experimental interest. In the context of lipid rafts |i)9^, confocal fluorescence 
microscopy studies have revealed the coexistence of fluid and solid domains on gi- 
ant unilamellar vesicles made of lipid mixtures. The shapes of these solid domains 
include stripes of various widths and orientations [100| IIOH I102| 1103] . 



3.2. Crystals of point particles 

Consider a collection of classical point particles constrained to a frozen (non- 
dynamical) two-dimensional surface S embedded in three-dimensional Euclidean 
space. The particles interact through a general potential defined in the three di- 
mensional embedding space or solely within the 2D curved surface itself. We focus 
primarily on the potential 

V{r) = . (58) 

Here, e is an "electric charge" such that if r is some quantity with dimensions of 
length then jr'^ has dimension of energy. The case 7 = 1 corresponds to the 
Coulomb potential in three dimensions. Although we do not treat this problem 
explicitly here, the replacement 

y(r) ^ -(r-^-1) (59) 
7 

allows us to treat the two dimensional Coulomb potential by taking the limit 7 — > 0, 

V{r) -e^\Qg{r) . (60) 
The electrostatic energy of a system of M particles at positions r{€)^ interacting 



via Eq.(58), with £ = (Zi, /2), ^1, ^2 G Z, becomes 



M 2 

Note that with this definition the power law interaction acts across a chord of the 
sphere, as would be the case for electron bubbles in helium. Although we focus here 
on curved crystals, there are some quantities which are insensitive to the curvature 
of the surface, and the simpler geometry of the plane can be used to compute them. 
The following two subsections hence focus on planar crystalline arrays of particles 



interacting via the potential Eq.(58). 
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3.2.1. Planar Crystals 



The electrostatic energy Eq. 



(61) and the corresponding elastic tensor, from 
may be explicitly computed 



which follows the elastic constants of the system 
for crystalline orderings of particles in a triangular lattice. For any non-compact 



surface, like the plane, the energy Eq. (61) is divergent for all 7 < 2. If 7 > 0, 
the divergence comes exclusively from the zero mode G = associated with the 
thermodynamic limit of infinite system size. This term (which would be subtracted 
off if a uniform background charge were present) can be isolated by setting \G\ = 
e <C 1 for this mode. The final result for the ground state energy reads 



2Eo 



+ Me 



r(7/2) {ac) 
2 ^ r(i-i) 



7(2 - 7) 
22-7 



a(7) 



lim 
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7/2 



Ac J 



vr 



r(i 



A 



c 



lim 



22-7 



7/2 \G\~.e • 



(62) 



and r 



Ac is the area of the unit cell of the triangular Bravais lattice {Ac — 2 
is the Euler Gamma function. The coefficient a is a sum over Misra functions. The 
coefficient ^(7) parameterizes the nonsingular part of the energy; its dependence 
on the exponent 7 is shown in Table [T} This negative quantity parameterizes the 
binding energy of the triangular lattice after the positive "zero mode" contribution 
is subtracted off. For 7 = 1, we have a two dimensional "jellium" model. In the 
case of a sphere no neutralizing background is necessary and the energy is rendered 
finite by the compact nature of the sphere itself. The maximum distance between 
points on the sphere provides a natural infrared cut-off. 
For small displacements of the particles from their equilibrium positions, one has 



E- 



1 



1 



\r{l) + u{t) - r{i') - u{£')\^ \r{i) - r{£')\^ 



(63) 



where u(i) is a small displacement of the particle £ in the plane of the surface from 
its equilibrium position r{i), and therefore a tangent vector to the surface. The 



elastic tensor Iljj (£,£') is defined as the leading term in an expansion of Eq.(63) 



(64) 



In deriving Eq.(63), we assume a constraint of fixed area per particle, enforced by 
a uniform background charge density or boundary conditions. This eliminates the 
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Tabic 1. Coefficients of the response func- 
tion Eq. ( |65[ l and the energy Eq. l |62^ . Re- 
sults are accurate up to six digits. The co- 
efficient p is a rational function of 7. 
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term linear in Ui{i). The final result is 

_ 22-Tvrr(l-7/2) piPj 
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r(7/2) 



12-7 
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7/2 
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+ 0{\p\' 



(65) 



The coefficients 7/(7) and p{j) depend only on the potential. In Table[T| some values 
of the coefficients for a range of potentials with < 7 < 2 are listed. 

3.2.2. Continuum free energy 

When the deviations from the ground state are small, the long wavelength lattice 
deformations may be described by a continuous Landau elastic energy of the form 



( 35 1 . The elastic tensor Eq. ( 64 ) , within Landau elastic theory, is then 
e^Iiij{p) = Ac [p\p?'5ij + (A + p)piPj\ . 



(66) 



A comparison with Eq. ( 65 ) immediately yields an explicit expression for the elastic 
constants of the crystal 



p = r/(7)- 



1+7/2 



00 



Y 



Ac 

4/^(A + p) 
2fi + X 



47?(7)- 



A 



1+7/2 ' 
c 



(67) 
(68) 



where Y is Young's modulus. The result A = 00 is equivalent to a divergent com- 
pressional sound velocity as p — > and for 7 = 1 is just a statement of the 
incompressibility of a two-dimensional Wigner crystal. Alternatively, we can allow 
for wavevector-dependent elastic constants p{p) and A(p) in Eq. (66). In this case 
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Figure 10. (Color online) Construction of an (n, m) icosadeltahedral lattice. The filled circles indicate two 
nearest-neighbor five fold disclinations. Because these defects sit on the vertices of an icosahedron, they 
are separated by a geodesic distance JJcos~^(l/-v/5), where R is the sphere radius. 



A(p) diverges as p ^ 0: 



Kp) 



22-')'7rr(l -7/2) 1 
r(7/2) ^ 



(69) 



while limp_^o l^{p) is given by Eq. (67). 

3.2.3. Spherical Crystals 

Spherical crystals have many properties not shared by planar ones, one of the 
most remarkable being that there is an excess of twelve positive (5— fold) disclina- 
tions. These disclinations repel, and the simplest spherical crystals will be those 
having the minimum number of defects (i.e. twelve) located at the vertices of an 
icosahedron. Triangular lattices on the sphere with an icosahedral defect pattern 
are classified by a pair of integers (n, m), as illustrated in Fig. 10 The path from 



one disclination to a neighboring disclination for an (n, m) icosadeltahedral lattice 
consists of n straight steps, a subsequent 60° turn, and m final straight steps. The 
geodesic distance between nearest-neighbor disclinations on a sphere of radius R is 
d = Rcos~^{\/ ^/b). The total number of particles M on the sphere described by 
this (n, m)-lattice is M = Nnm, where 



10(n^ + + nm) + 2 



(70) 



Such (n, m) configurations were once believed to be ground states for relatively 
small numbers (M < 300, say) of particles interacting through a Coulomb potential 
|104| I105| 11061 1107| 1108] . The energy of discrete particle arrays described by Eq. 
(61) can be evaluated by starting with some configuration close to an (n, m) one 



and relaxing it to find a minimum. It is found that the (n, m) configurations are 
always local minima. Whether these icosahedral configurations are global minima 
as well will be analyzed later. Results for the energy E{M) are shown in the inset 



to Fig. 11 



From Fig. [IT] it is clear that energies grow very fast for increasing volume. More 
interestingly, the (n, 0) and (n, n) configurations show a growing difference in en- 
ergy for increasing volume size, implying that the energy of icosahedral configura- 
tions does not tend to a universal value for large numbers of particles but rather 
remains sensitive to the (n, m) configuration, a result also noted by other authors 
|105] . Further insight comes from investigating the distribution of energy. Plots of 
the local electrostatic energy, the electrostatic energy at point x on the sphere, as 



defined in Eq. (61 ) are shown in Fig. 12 |109j . From Fig. 12 it should be noted that 
the triangles obtained by the Voronoi-Delaunay construction, after minimization 



of the potential Eq. (61 ), are very close to equilateral. The distribution of the local 
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Figure 12. (Color online) Potential energy distribution for a (ri, 0) configuration with n=10 and M=1002 
(top) and a (n,n) configuration with n=6 and M=1082 (bottom). 

energies for the (n,0) and {n,n) configurations are very different. The (n,0) con- 
figuration shows maximum energies along the paths joining the defects. The (n, n) 
configuration, on the other hand, has its maximum energies along the directions 
defined by the triangles formed by three nearest neighbor disclinations. The size of 
these regions of differing electrostatic energy turns out to scale with system size, 
making it very plausible that there might be small differences in the energy per 
particle for (n, 0) and (n, n) configurations in the limit n oo. 



3.3. Geometric Formalism on the Sphere 

Spherical substrates provide the simplest example of the problem of crystals on 
curved surfaces. The study of spherical crystals is simplified by two important 
properties: there is a unique scale with dimensions of length, the radius R, and 
there is a fixed excess disclinicity of twelve following from the Euler theorem ([T]) 

N 

Q = ^% = 6x = 12. (71) 

i=l 



The free energy Eq. (55 ), applied to the sphere, is tractable analytically because the 
inverse biharmonic operator on a sphere of radius R can be computed explicitly. It is 
shown in [73] that the Green's function for the biharmonic, in spherical coordinates 
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has the fohowing simple form on a unit sphere: 



1-C05/3 

x(r,(/.'^;0^/) = l+ / ' , (72) 

JO I — z 

where /3 is the geodesic distance between two dischnations located at {9"',(t)"') and 

(Qb Ab\ 



cos /3 = cos d'' cos 9'' + sin 0" sin cos((/>'' - /) . (73) 
The total energy of a spherical crystal with an arbitrary number of dischnations 



follows from Eq. (55) and Eq. (72) and has the simple form [73 



2E{Y) = Eo + —R" mx{0\ 0^ e^,<j>^) + (74) 

i=l j=l 

where {6i, (pi}i=i^... are the coordinates of defects and we restrict ourselves to 
5— fold {qi = +l)and 7— fold (qi = —1) defects. The quantity Eq is the zero point 
energy and is defined in Eq. ( |62[ ) . Although 5 and 7- fold dischnations will in general 
have different core energies [llOj , we assume equal core energies here for simplicity. 
What matters for our calculations in any case is the dislocation core energy E^, 
which we take to be E^ = E^ + E-j = 2Ec. 

The value of the Young's modulus and the flat space ground state energy Eq 



have been computed in Sec. 3.2.1 When the sphere radius R is large compared to 
the particle spacing a, we can use flat space values of Y and the flat space energy 
Eq{M) associated with a finite number of particles M. To obtain the leading terms 
in the expansion of the ground state energy for large but finite M, the precise 
compactification of the plane employed is irrelevant - it may be achieved by periodic 
boundary conditions, for example. For a sufficiently large plane the finite size effects 
will be negligible. The density a of particles is then M divided by the total surface 
of the compact plane, taken to be the surface area 5 of a sphere of radius R, 

a = l/Ac = ^ ^ S = AnR^ . (75) 



From Eq. (67) the expression for the Young's modulus suitable for M particles on 



a spherical crystal of radius R with < 7 < 2 is then 



4,(7)M^+^/^ ^ (76) 
^ (47r) 1+7/2 R2+y ' ^''^> 



One remaining detail is the divergent contribution to the energy £^0 in Eq. (62). 
Since the divergent part comes solely from the zero mode, the spatial variations in 
the density of the actual distribution are irrelevant. It may therefore be computed 
for a uniform density of charges. The divergent part is identical to the energy of 



a constant continuum of charges as described by the density Eq. (75). We now 
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evaluate this divergent part of the energy on a sphere, instead of a plane. 
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The divergent part has thus been regularized, and the energy is finite and well- 
defined for all M < 00. Note that for the case -/ < 2, Ed M'^-^/'^{M/ S)^'"^ . 
Hence Ed is not simply a function of the particle density M/S, as one would 
expect for a short range interaction. 

3.3.1. The energy of spherical crystals 



Upon substituting the elastic constant of Eq. (76) into Eq. (74), one arrives at 

Try ^ ^ 



i=i j=i 
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+ NE, 



(78) 



where Eq is defined in Eq. (62 ) and the function C{ii ■ ■ ■ in) depends on the position 
H = (^ii 01 ) etc. of the N disclination charges and is universal with respect to the 
potential. The total energy of a spherical crystal, including the contributions to Eq 
is then 
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Note that the leading correction to the zero mode energy proportional to 
varies as M^"'"'''/^, and depends both on the flat space function ^(7) and on the 
C-coefficient 



N N 

C(ii,-- - ,iiv) = EE'?*«^-^(^''^';^''^') ' 
j=i i=i 



(80) 



associated with a particular configuration of disclinations. Note that the core en- 
ergies contribute to the second sub-leading coefficient. For short-range potentials, 
such as 7 > 2, the ground energy is extensive, and the leading term varies as 
7^^1+7/2 _ ^iiQ extensive nature of the M^'^^^'^ term becomes clear upon noting that 



R-r 



(81) 



where a is the particle spacing. Comparison with Eq. (74 ) shows that the dimension 
of Young's modulus Y arises solely from the lattice constant a and the electric 
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Figure 13. (Color online) Illustration of the calculation discussed in the text. The energy of the discrete 
(n, n) configuration on the left is extrapolated for large M and compared to the energy computed with 
the continuum model on the right. While in the continuum model only twelve degrees of freedom (the 
12 disclinations) need to be considered, the direct calculation of a family of discrete models requires the 
consideration of the full lattice and a careful extrapolation of the energies to large M. 



charge e, consistent with elastic constants arising from physics on the scale of the 
lattice constant in an essentially flat geometry. This observation is now generalized 
to the rest of the couplings discussed in the previous section. 

For a fluid membrane not on a frozen topography, Helfrich terms arising from 
the mean curvature H as well as the Gaussian curvature can be important. Their 
scaling behavior is given by 



Kg 



j d^xK ^ K j d^xH^ ^ ^M^/"^ . (82) 



Both terms would therefore contribute to the same order in the M expansion, 
although the last term is purely topological. For crystals embedded in a frozen 



topography we expect an expansion along the lines of Eq. ( 79 ) 



e2 



2EUM) = \ aoM^ - ^ a.M^/^+i-^ ) — . (83) 
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The nonextensive term a^M"^ arises from the long range interactions. The next 
extensive contribution comes from the interaction between Gaussian curvature and 
defects as well as the extensive energy per particle in flat space. Bending rigidity 
contributions are higher order in 1/M and can be absorbed into a redefinition of 
the disclination core energy. Core energies also depend on non-universal details of 
the short-distance physics. 

The results presented so far are strictly for systems at zero temperature. 
In systems with short range interactions, the elastic constants can be strongly 
temperature-dependent. An extreme example is hard disks of radius ao, which may 
be viewed as a limiting case of a power law potential of the form V{r) ~ eo(ao/?')''', 
with 7 — > oo. In this case, the elastic constants are strictly proportional to tem- 
perature. It is straightforward, however, to adapt the techniques presented here to 
the simpler problem of short range pair potentials. 

3.3.2. Energies of Icosahedral configurations 

The configuration on a sphere with the minimum number of charge ±1 defects is 
twelve -|-1 (5— valent) disclinations, which minimize their energy by sitting at the 
vertices of an icosahedron y. The energies of such configurations will be computed 



for the discrete spherical tessellations described in Sec. 3.2.3 and compared with 
the predictions of continuum elastic theory, as illustrated in Fig. [TS} It is well 
established that for sufficiently large values of M configurations with more than 12 
disclinations (i.e., those with "grain boundary scars") have lower energies |73 |llll] . 
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Figure 14. (Color online) Energy coefficient ai as a function of gamma (solid line) and from the numerical 
results with (n,m) configurations (ffiled circles), for the icosahedral configurations. Plot of ai(7,y) - 
ai(7,M^"'"^ (circles), ai{'y,y) - ai(7,y)(">0) (diamonds), ai (7, 3^)("'") - ai(7,y)("'0) (squares). 



It is of interest, however, to study simple icosahedral configurations for large M, 
as metastable states with a well defined energy. 

Within the continuum elastic theory it can be shown that twelve disclinations 
at the vertices of an icosahedron minimize the energy [73] when no further defects 



are allowed. The C-coefficient of Eq. ( 79 ) for this configuration of defects has been 
computed in [TH] 



C{y) = 0.6043 . 



(84) 



y here stands for a particle configuration with 12 defects at the vertices of an 
icosahedron. Using the energy of Eq 
the expansion of Eq 



(79) 



the coefficient 01(7,3^) appearing in 
(83) may be computed, with the results shown in Fig. 14 



and Table [2] From the results described in Sec. 3.2.3 the ai coefficient may be 
extrapolated to very large numbers of particles using the expansion derived from 
Eq. (83). Indeed, as shown in Fig. 15 plots of 



6(M)^ 



(85) 



vs 1/M are linear, with a slope that determines 01(7) and an intercept related 
to the higher order core energy-like contribution. The results of these extrapola- 
tions are shown in Table [2} The agreement between the continuum elastic theory 
and the explicit computation for the (n, n) configuration is remarkable, holding 
to almost five significant figures. For the (n, 0) lattice there is agreement to four 
significant figures. This agreement is even more striking when it is recalled that the 
ai-coefficient is obtained after subtraction of the term ao{'y)M'^, as illustrated in 
Fig. 11 Furthermore, in the range from 7 = 0.125 to 7 = 1.875, all the significant 
digits vary and yet the accuracy of the calculation is virtually independent of 7. 

3.3.3. The Energy difference of the {n,m) lattices 

The oi-coefficient computed within the continuum elastic approach above does 
not depend on the icosadeltahedral class (n, m). Results from the direct minimiza- 
tion of particles do, however, show a weak dependence (in the Ath significant digit) 
on the particular {n,m) configuration, as is apparent from Fig. 11 and Table |2} It 
should be noted that the discrepancy from the continuum result has a well defined 
sign, and is therefore reasonably attributed to a term not present in the energy 
expansion. 
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Figure 15. (Color online) Numerical estimate of e(AI) as a function of 1/M for (n, 0) and (n, n) icosadelta- 
hedral lattices with 7 = (1.5,0.5). 



Table 2 . Numerical values of the coefficient 
^11(753^) (twelve disclinations on the vertices of an 
icosahedron) using the C-coefficient from Eq.(|S4j. 
The same coefficients from the (n, n) and (n, 0) 
lattices 



7 




in, n) 


(n,0) 


1.875 


4.45118 


4.45110(4) 


4.45095(4) 


1.75 


2.47175 


2.47166(3) 


2.47150(3) 


1.625 


1.82629 


1.82621(2) 


1.82603(2) 


1.5 


1.51473 


1.51454(2) 


1.51445(2) 


1.375 


1.33695 


1.33683(4) 


1.33667(4) 


1.25 


1.22617 


1.22599(7) 


1.22589(7) 


1.125 


1.15366 


1.1535(2) 


1.15340(2) 


1.0 


1.10494 


1.10482(3) 


1.10464(3) 


0.875 


1.07187 


1.07174(3) 


1.07160(3) 


0.75 


1.04940 


1.04921(6) 


1.04910(6) 


0.625 


1.03421 


1.03413(5) 


1.03398(5) 


0.5 


1.02392 


1.02390(4) 


1.02372(4) 


0.375 


1.01672 


1.01663(6) 


1.01656(6) 


0.25 


1.01115 


1.01106(3) 


1.01103(3) 


0.125 


1.00595 


1.00592(2) 


1.00589(2) 



3.4. Thomson problem with a continuous distribution of dislocations 

When the number of particles is extremely large, the minimum energy configura- 
tions can be approximated by a closed analytical form, upon assuming a continuous 
distribution of defects. 

The formal elimination of the geometric frustration introduced by the Gaussian 
curvature may be formulated as a concrete set of equations in the case of the sphere. 
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We shall use the identity 



vr 



N 



i=l 



oo I 



N 



(86) 



oo I N 
1=1 m=—l i=l 



1=1 m=—l 



i=l 



which follows from the topological constraint Eq. (71). Provided a disclination 
configuration exists such that 



AT 

E«^^- 
1=1 



i) = 



(87) 



for each (/ > 1,iti), the disclination density completely screens the Gaussian cur- 



vature. A configuration of defects satisfying Eq.(87) is an absolute minimum of the 



elastic energy, a result easily understood by writing the energy in the form 



E 



^2y 
Eo + —R 



oo 

1=1 m=—l 



\l^i=l 'ii^r 



+ NE, 



where the zero point energy Eq in Eq. (62) is kept. A configuration satisfying Eq. 



(87)) will be denoted by Q. For this hypothetical configuration, the C-coefficient 



in Eq. 79 vanishes, although there is now a large contribution (linear in R) from 
the dislocation core energies represented by the last term of Eq. ( 88 ) . 



Table 3. 


Value of 


the ai coc 


fficients for 


the Q configuration 


Eq. |87]l. 




7 


ai{7.G) 


1.875 


4.45227 


1.75 


2.47289 


1.625 


1.82746 


1.5 


1.51592 


1.375 


1.33815 


1.25 


1.22737 


1.125 


1.15485 


1 


1.10610 


0.875 


1.07297 


0.75 


1.05044 


0.625 


1.03515 


0.5 


1.02473 


0.375 


1.01737 


0.25 


1.01161 


0.125 


1.00620 



The Q configuration may be characterized more explicitly. It consists of a density 
of dislocations that converges to the local Gaussian curvature. It can be shown 
that upon approximating the dislocations (each regarded as a disclination dipole 
with spacing a) as a continuum distribution, this dislocation density for a sphere 
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becomes 



fc=i 

The summation here runs over the six coordinates of the northern hemisphere of 
an icosahedron, (0,0) and {6y, 27rk/5), where 9y = arccos(l/\/5) and is the 
angle 9 relative to a coordinate system with the north pole located at {9y, 27rk/5) 
for k = 1, ■ ■ ■ 5. This angle is given implicitly by 

cos[afc(6', 4>)] = cos{9) cos{9y) - sm{9) sin(6'3;) cos (^k + . (90) 

The implicit form of Eq. ( 89 ) can be further simplified 

= f\9, 0) I - sm{9y) sin (^^k + 0^ ge 

+ [cos{9y) sin 9 + sm{9y)] cos 9 cos ^0 + 



(91) 



where f'^{9^(f)) = l/sin[a,fc(0, (/>)]. Close to one of the 12 disclinations with charge 

(92) 



s = 2tt/Q Eq. (89) predicts a singularity in the dislocation density ^74j 

s 



2'KRa 

For small angles, close to each disclination, there is a short-distance singularity 



in agreement with known results in flat space. Eq. (89) represents a continuous 
distribution of dislocations, and neglects both dislocation discreteness and their 
mutual interactions. It represents six families of dislocations with azimuthal Burg- 
ers vectors associated with antipodal pairs of the 12 original disclinations in the 
icosahedron. When discreteness and interactions are taken in account, we expect 
these dislocations to condense into grain boundary arms, containing with quan- 
tized Burgers vectors and variable spacing in the radial direction \7?>\ \\\2\ 1113] . 
The total number of discrete arms remains, therefore, the variable that needs to 
be determined for a discrete solution of the Thomson problem. 

3. 4-1- The intermediate regime 

Within the continuum elastic approach, the dominant configurations for a small 
number of particles are 12 defects with an icosahedral symmetry [73]. We have 
just seen, however, that adding a continuous distribution of dislocations, as might 
be appropriate when the particle number is large, can more efficiently screen the 
Gaussian curvature on a sphere. The natural problem then becomes to determine 
the precise structure of the defect arrays for intermediate numbers of particles 
when the discreteness of interacting dislocations is taken into account. 

We note first that the particular arrangement of defects dominating in this regime 
will not be fully universal. The particular array structure favored can vary from 
system to system with fixed particle number, depending, e.g., on details such as 
the dislocation core energy. This result may be traced back to the M-expansion of 
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Figure 16. (Color online) Results of a minimization of 500 particles interacting with a Coulomb potential, 
showing the appearance of scars. 



r • 
■ 



Figure 17. (Color online) Ground state configurations for M 2000 particles obtained from the continuum 
elastic formalism. In the top figure one finds scars (m = 2) and in the bottom pentagonal buttons (m = 5) 
forming a rhombic tricontahedron. 



Eq. ( 83 ) , in which the sub-leading terms which depend on non-universal properties 



influence the dominant terms for finite values of M. Some typical defect configura- 
tions obtained by direct minimization of particles on the sphere are shown in Fig. 
16 and show incipient scars, already at 500 particles (in [73] the minimum number 
of particles where scars are systematically found is predicted to be around 400). By 
using the geometrical model described here, where the energy is parameterized just 
by a Young's modulus and a dislocation core energy |73| |92] one can simulate larger 



particle numbers and one obtains results as in Fig. 17 Note the occurrence of low 
energy configurations with scars (m = 2) in one instance and pentagonal buttons 
(m = 5) in another. The dislocation spacing decreases the further a dislocation is 
from the central disclination. 



An overview of results involving grain boundary scars is presented in Fig. 18 



If a disclination is placed on a perfect crystal, no additional defects will appear 
if the disclination is located on the tip of a cone with total Gaussian curvature 
equal to the disclination charge. If a disclination is forced into a flat monolayer. 
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Figure 18. (Color online) Schematic illustrating the genesis of grain boundary scars. A disclination is first 
constructed from a perfect lattice. If this disclination is placed on a tip of a cone, with a delta function of 
Gaussian curvature balancing the defect charge, then no additional defects form. If the crystal is forced 
into a monolayer, grain boundaries radiating out of the disclination radiate all the way to the boundary. 
In the intermediate regime of constant non-zero Gaussian curvature, m grain boundaries of finite length 
and variable spacing of dislocations form. 



then m low angle grain boundaries, with constant spacing between dislocations as 
shown in Fig. 18 and grains going all the way to the boundary, will be favored 
(see |112j for a detailed discussion). In the intermediate situation where a finite 
Gaussian curvature is spread over a finite area, as in the case of a spherical cap, 
a disclination arises at the center of the cap and finite length grain boundaries 
stretched out over an area of (7r/3)i?^ with variable spacing dominate, again as 



illustrated in Fig. 18 



Additional results may be obtained for the number of arms within the grain 
boundary, the actual variable spacing between dislocations within the grain and 
the length of the grains as a function of the number of particles. 

When grain boundary scars appear, one can estimate the number of excess dislo- 
cations which decorate each of the 12 curvature-induced disclinations on the sphere 
using ideas from Ref. ^3]. This estimate is in reasonable agreement with experi- 
ments probing equilibrated assemblies of polystyrene beads on water droplets [22] • 
Consider the region surrounding one of the 12 excess disclinations, with charge 
s = 27r/6, centered on the north pole. As discussed in Ref. [73], one expects the 
stresses and strains at a fixed geodesic distance r from the pole on a sphere of 
radius R to be controlled by an effective disclination charge 



^eff 



(r) 



JO 



vr 



47rsin^ f — ) • 
3 \2R) 



(94) 



Here the Gaussian curvature Is K = 1/E? and the metric tensor associated 
with spherical polar coordinates {r,(f)), with distance element ds^ = d'^r -\- 
s\v?{r / R)d^4), gives ^/g = i?sin(r/i?). Suppose m grain boundaries radiate from 
the disclination at the north pole. Then, in an approximation which neglects inter- 
actions between the individual arms, the spacing between the dislocations in these 
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grains is ^73] 



l{r) 



am 



which imphes an effective dislocation density 
11" 



nd{r) 



^ A ■ 2 ( 

47r sm — - 

3 \2R 



l{r) ma 

This density vanishes when r ^ rc, where 



27r 
ma 



cos 



(95) 



(96) 



Tc = i?arccos ( - ) Ri 33.56°i? , 



(97) 



which is the distance at which the m grain boundaries terminate. The total number 
of dislocations residing within this radius is thus 



Nd = m dr nd{r) 
Jo 



IT 

7 

3a 

vr 
3 



4tt 



11 — 5 arccos 



0.408 ( - 

a 



r 



dr sin^ ( „ 
a ./o ^2i? 



(98) 



As discussed in Ref. |114| . it is also of interest to consider 2tt disclination defects 
(appropriate to crystals of tilted molecules |115j ) on the sphere. The icosahedral 
configuration of 12 s = 2tt/6 disclinations is now replaced by just two s = 2it discli- 
nations at the north and south poles. Using the approximation discussed above, it 
is straightforward to show that the density of dislocations in each of m (noninter- 
actmg) grain boundary arms now reads 



nd{r) 



2tt / r 

cos - 

ma \K 



(99) 



This density vanishes at Tc = (7r/2)i?, corresponding to a hemisphere of area on 
the sphere for each cluster of arms. 

It is of considerable interest to repeat the above calculation for a square lattice, as 
found for example in the protein surface layers (s- layers) of some bacteria |116tlll7] . 
In this case the basic disclination has s = 27r/4. The effective dislocation density 
becomes 



nd{r) 
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(100) 



This density vanishes when r ^ re, where 



R arccos 



41.4°i? , 



(101) 
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which is the distance at which the m grain boundaries terminate. The longer angu- 
lar length of square lattice scars reflects the larger initial disclination charge (90°) 
that must be screened. The total number of dislocations residing within this radius 
is thus 



Thus the angular length of scars and the total number of excess dislocations is a 
measure of the underlying topology of the lattice tiling the sphere. 

3.5. Interstitial fractionalization on the sphere 

3.5.1. Introduction 

Consistent with theoretical expectations, experiments on particle-coated water 
droplets in oil [92] reveal that the twelve excess disclinations sprout grain bound- 
ary "scars" for sufficiently large R/a. When triangulations of microscopic particle 
packings are used to reveal the local coordination number, these grain boundaries 
appear as additional dislocations, i.e., 5-7 pairs, arrayed around an unpaired 5, in a 
pattern such as 5-7-5-7-5-7-5-7-5. Although the critical value of R/a above which 
grain boundaries appear depends on microscopic details, both theoretical estimates 
and experiments indicate that this instability arises as soon as R/a ^5 — 6, i.e., 
when the total number of particles exceeds several hundred. Thus, unlike crystals 
in flat space, dislocations arrayed in grain boundaries are an intrinsic part of the 
ground state. These grain boundaries can, moreover, stop and start freely on the 
sphere, unlike their flat space counterparts. Such terminations occur naturally (and 
with low energy cost) because crystalline grains rotate under parallel transport due 
to the nonzero Gaussian curvature of the sphere. 

If disclinations and dislocations are crucial for understanding spherical crystal- 
lography, what can we say about vacancies and inter stitials, which are well known 
to play a key role in conventional crystals [1181 1119] . It is natural to introduce 
vacancies and interstitials in an attempt to understand spherical particle packings 
that deviate from the sequence of "magic numbers" Nnm = W{n? + nm + m?) + 2. 
It is tempting to ignore the instability to grain boundary scars for large 
and regard the commensurate (n, m) lattice as an interesting metastable state. It 
would then be natural to introduce vacancies and interstitials to describe candi- 
date ground state packings for Nnm i t particles, where t = 1,2, ... and much less 
than the distance to the next magic number. There is, however, another surprise 
in store: in contrast to flat space, where vacancy and interstitial defects are sta- 
ble and well defined, one finds that interstitials and vacancies are typically ripped 
apart into dislocation fragments by the strain fields of nearby 5-fold disclinations. 
The dislocations then combine with some of the excess 5's to form defect clusters 
such 5-7-5. Thus, vacancies and interstitials lose their integrity via fragmentation 
in spherical crystals, and mediate formation of small grain boundary scars. 







(102) 
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Figure 19. (Color online) Initial icosahedral configuration of an (8,3) spherical crystal lattice with 12 
disclination defects. This lattice is chiral, and is arranged such that 8 steps along a Bragg row, followed by 
3 steps to the right along a Bragg row at a 120° angle, connect neighboring 5-fold disclinations. Distinct 
initial locations for interstitials are shown as triangular plaquettes labeled by characters and numbers. 



3.5.2. Interstitial fractionalization 

We are interested in studying the effect of inserting an interstitial or a vacancy 
into a regular (n, m) icosahedral lattice by adding or removing a single particle, 
giving rise to particle numbers A^^^il not falling in the classification of icosadelta- 
hedral lattices and in studying the relation of interstitials to the extra dislocation 
defects (scars) found in spherical crystals above a critical system size [73l [92] . 
Refs. [I2Q1 [I2ll [ISa una UMI ESSI EZT] study configurations and energies of 
interstitial and vacancy defects and their energetics in triangular lattices in flat 
space. 

The presence of excess dislocation defects in the ground state of spherical crystals 
is dramatically illustrated by the following numerical experiment: we start with a 
regular icosadeltahedral tessellation of the sphere - say an (8, 3), corresponding to 



Ns3 = 972 (Fig. 19). This may be done with the applet located at |109j using the 
Construct (m, n) algorithm. Although the true ground state for 972 particles on 
the sphere with most pair potentials has additional dislocation defects (i.e. tightly 
bound pairs of 5- and 7- coordinated particles) arrayed in grain boundary scars [73] , 
the regular icosadeltahedral lattice is a local minimum from which it is difficult to 
escape without the addition of thermal noise. In fact it is a major challenge to find 
fast and reliable algorithms to locate the true ground state (global minimum) in this 
problem with its complex energy landscape. Now add a single particle to the lattice 
at the center of mass of a spherical triangle whose vertices are 3 nearest-neighbor 
5-fold disclinations (using commands shift -|- click). The self-interstitial so formed 
is then relaxed by a standard relaxation algorithm, with sufficient thermal noise to 
allow dislocation glide over the Peierls potential |93] . One immediately finds that an 
interstitial is structurally unstable. In a few time steps it morphs into a complex of 
dislocations with zero net Burgers vector. The most common structure observed is 
a set of three dislocations, with Burgers vectors perpendicular to a line joining each 
5-7 pair, inclined at 120° angles to each other. Eventually the interstitial complex 



is ripped apart entirely, as illustrated schematically in Fig. 20 (see also Fig. 21). 
Intermediate configurations and final states as the dislocations glide apart will be 
classified later. Most often three separate dislocations are formed which each glide 
toward a 5-fold disclination. The end result is the formation of a "mini-scar" (a 
5 — 7 — 5 grain boundary) at each of the vertex 5s. Subsequent removal of a particle 
to restore the particle number to the original 972 and relaxing still leaves scars 
with total energy lower than the starting configuration with 12 isolated 5's. This 
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(a) (b) (c) (d) 

Figure 20. (Color online) A schematic of interstitial fractionization. The H symbols are alternative ways 
of representing dislocations depicted elsewhere as 5-7 pairs. 



observation confirms that scars are definitely lower energy states and not simply 
artifacts of the relaxation algorithm. 

The above phenomenon of low-temperature (T > 0) unbinding of dislocations by 
spatial curvature is a curved space analog of melting at finite temperature. The 
extended nature of fractionated interstitials (each separating dislocation compo- 
nent involves an extra row of particles) means that they cannot be treated as small 
perturbations from the initial spherical crystal with particle number Nnm- 

Let's return to the specific case of the = 973 particle configuration generated 
by an interstitial inserted into one triangular plaquette of a regular (8, 3) lattice 
of Nnm = 972 particles with the requisite 12 disclination defects (5-fold coordi- 
nated particles) at the vertices of a regular icosahedron, as shown in Fig. 19 The 
spherical crystal is distorted by the additional particle - the local configuration 
adopted by the interstitial changes as the crystal relaxes toward a lower energy 
state. As in the case of planar lattices, one also finds here that various intersti- 
tial defect configurations appear, such as the twofold symmetric interstitial I2, the 
threefold symmetric interstitial Is, and the fourfold symmetric interstitial I4 (see 



Figs. 21, 22 and 23). The most common intermediate complex formed by the 



interstitial is threefold symmetric in the rough form of a triangular loop composed 
of three dislocations with radially oriented Burgers vectors. All of the configura- 
tions adopted by an interstitial prior to unbinding can be described as a set of 
dislocations with zero net Burgers vector. 

In marked contrast to interstitial defects in a planar crystal, interstitial defect 
configurations in a spherical crystal are metastable states with characteristic decay 
processes. As we shall see, the instability is caused by interactions with the in- 
evitable disclinations associated with the nonzero Gaussian curvature of the sphere. 
A representative evolution of an interstitial inserted at the center of three neigh- 
boring disclinations in an (8, 3) icosadeltahedron is shown in Fig. 21 After some 
local relaxation the interstitial configuration denoted I3 is formed, as shown in 
Fig. 21 (a). We also show there the construction of a triangular contour surround- 
ing the original interstitial defect. The presence of an interstitial follows because 
the contour encloses 7 particles instead of 6, the number appropriate to a per- 
fect triangular lattice. During the annealed relaxation illustrated in Figs. 21 'b) 
through (d), the dislocations which are bound together in the initial interstitial 
unbind into individual dislocations and subsequently glide towards nearby discli- 
nations, eventually forming minimal 5-7-5 grain boundary scars. If one thinks of 
the initial dislocations as internal degrees of freedom within the interstitial, one 
could say that one-third of an interstitial is present in each mini-scar in the final 
state and thus that the interstitial demonstrates 1/3 fractionalization. For other 
initial conditions the fractionation is into two dislocations, each representing 1/2 
of the original interstitial. The instability of interstitial defects in curved crystals 
may be studied via continuum elasticity theory by calculating the interaction en- 
ergies between defects at each stage of the relaxation process of Fig. 21 |119] . We 
note that the triangular plaquette around the initial defect has been deformed such 
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Figure 21. (Color online) The fractionalization of an interstitial defect inserted at the center of three 
neighboring disclinations in an (8,3) spherical tessellation: (a) the initial interstitial defect configuration 
(73) shown surrounded by a triangular reference contour; (b) the bound defect unbinds to form three 
separate dislocations (5-7)- pairs; (c) each dislocation glides (i.e., moves parallel to its Burgers vector) 
towards the nearest disclination; (d) three mini-grain boundary scars are formed. We also keep track of 
the evolution of the initial defect in (a) by illustrating the deformation of the triangular contour around 
the initial defect for (b), (c), and (d) induced by the passage of the dislocation. 

that it conforms as closely as possible to the regular triangular lattice during the 
relaxation process. Its deformation reveals the passage of the escaping dislocations. 

3.5.3. Interstitial defect energies 

In this section, we discuss the energy of an interstitial defect in a spherical crystal. 
Consider point particles constrained to lie on the two-dimensional surface of a 
unit sphere. The energy of N particles interacting through a generalized Coulomb 
potential within the curved surface is given by 



E 



l,N 



1 



(103) 



where x is the position of the particle in three dimensions and 7 is an integer. For 
a flat triangular lattice with periodic boundary conditions, the interstitial defect 
energy at constant density was defined in |122| 1123] as 



El = E, 



relaxed 



E. 



per ) 



(104) 



where E^eiaxed is the relaxed energy of the rearranged lattice of N particles with 
the interstitial defect in the area A, and Ep^r is the energy of the perfect crystal 
at the same areal density N/A. In curved space, however, the definition of a "per- 
fect crystal" is more subtle, since disclination defects resulting from the Gaussian 
curvature and the topology are inevitable. We will take as a reference crystal the 
(n, m) icosadeltahedral configurations corresponding to triangular tessellations of 



Two-Dimensional Matter: Order, Curvature and Defects 



43 



a magic number of particles Nnm = 10(n^ + m? + nm) + 2. Once an interstitial 
or vacancy is added to such a (n, m) configuration, we are no longer at a magic 
number of particles since these are quite sparsely distributed. We thus need to de- 
fine the energy of the perfect crystal. Here we define the energy of the interstitial 
(vacancy) defect at constant density in the spherical crystal as 

El = Elocal - E^nnealed i (105) 

where Eiocai measures the energy of the relaxed interstitial while the constituent 
dislocations are still bound and ^a„„eaied minimum energy of all possible 

final states attained after annealed relaxation leading to interstitial fractionaliza- 
tion. This definition will be more explicitly discussed in the following section (see 
Table [5| . We note that both Eiocai and -E'^^^eaied measured at the same areal 
density {Nnm^'^)/ A, where ±1 correspond to an interstitial (vacancy) respectively. 
The lowest relaxed energy E*^^^^^^^ plays the role of the energy of the perfect lattice 
in the planar case at the density of {Nnm i 1)/A. 

Numerical measurements of Eiocai and E*^^^^^^^ for the power- law potentials with 
7 = 1, 3, 6 and 12 have been performed by adding one interstitial at the center of a 
spherical triangle formed by three nearest-neighbor disclinations in the (8, 3) lattice 



(the location represented by C in Fig. 19) |119] . Eiocai is measured by quenching 



the system at the moment just prior to the fractionation of the interstitial into 
individual dislocations [(a) in Fig. 21 . The results are reported in Table |4j 



Table 4. The lowest local and annealed relaxed energy with the central interstitials 
created by putting a particle at C in Fig. |19| of the (8,3) lattice, for 7 — 1, 3, 6, and 
12. The differences between two relaxed energies are calculated. Because the particles 
are embedded in a sphere of unit radius with our conventions, near-neighbor particle 
spacings are of order a ~ N~^^'^ <C 1, leading to a strong 7 dependence in the total 
energy given by Eq. jl03[>. 



7 


1 


3 


6 




12 


Local 


456601.99 


2840600.7 


9.62182 


xlO** 


3.03015 xlO" 


Annealed 


456600.91 


2840025.5 


9.60570 


xlO* 


3.00313 xlO^-* 




1.08 


575.2 


0.01612 


xlO* 


0.02702 xlO" 



3.5.4- Position dependence of interstitial defect energies 



By adding a particle in different plaquettes within the large spherical triangle 
of Fig. 



19 



one can explore the dependence of the final state on the initial inter- 
stitial location. In contrast to the case for planar crystals, both the location and 
orientation of the interstitial defect relative to nearby disclinations influences the 
resultant configuration and its corresponding evolution, leading to distinct relaxed 
configurations. The insertion of a particle at the center of the large spherical trian- 
gle leads to an /3-type initial configuration, whereas adding the interstitial to the 
plaquette along an edge results in an /2-type initial configuration. 

During the relaxation process, one also finds that the dislocation complex rep- 
resenting an interstitial can rotate so as to reorient the Burgers vectors so that 
constituent dislocations can glide towards a nearby disclination and bind to it. 
This phenomenon is especially noticeable if one places one extra particle slightly 



off the absolute center C, such as at the locations C or in Fig. 19 



In Figs. 22 and Fig. 23 this phenomenon is illustrated explicitly. In Fig. |22[a), 
an interstitial initially placed at the position C in Fig. 19 morphs quickly to an I3 



configuration. In (a), however, the orientations of the dislocations within are not 
appropriate for dislocation glide to the surrounding disclinations since the 5 end 
of the dislocations point towards rather than away from these 5-fold disclinations. 
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Table 5. The three classes of final annealed state depending 
on the initial interstitial location. The relaxed energies are 
measured for the power law potential 7 — 6. Here Eannealed 
with 3 scars corresponds to -E^aTiTieaZed ■ 



Initial Location 


-^annealed 


Annealed State 


C, C, 7, 8, 10 


9.60570 xlO* 


3 scars (5-7-5) 


0, 2, 3, 4, 6, 9 


9.61011 xlO* 


2 scars (5-7-5) 


1, 5 


9.61062 xlO* 


2 scars (5-7-5-7-5) 



causing them to be repelled. Glide-induced fractionalization is therefore prohibited 
in this orientation. Remarkably, though, the entire complex of dislocations can 
change its orientation by a transition through an intermediate I2 configuration 



(shown in Fig. 22 'b)) and subsequently to a second 13 configuration [Fig. [22[c)]. 
The final I3 configuration is rotated by 60° with respect to the first I3 and can 



now fractionate analogously to an interstitial with initial position C in Fig. 19 



One also find rotational reorientation of an interstitial defect with an 14-type in- 
termediate state, as shown in Fig. 23 This particular relaxation process reveals an 
interesting feature of dislocation dynamics on a curved surface. The I3 configura- 
tion generated after the intermediate [see Fig. [23|^c)] now has one dislocation with 
its glide plane such that it can glide head-on into a vertex disclination. This discli- 
nation absorbs the dislocation but hops over one lattice spacing to accommodate 
the curved space Burgers vector. The other two dislocations end up bound in the 
form of mini-scars. In this case then the interstitial has fractionated into 2 rather 
than 3 parts and one say that there has been 1/2 fractionization of the intersti- 
tial. Absorption and emission of dislocations by 5-fold disclinations are somewhat 
analogous to absorption and emission of vacancies and interstitials by dislocations 
(allowing dislocations to climb), a phenomenon well-known in flat space. 

What is the dependence of the final state on the initial location of the interstitial? 
The distinct initial conditions shown in Fig. 19 lead to three different final annealed 
states, as summarized in Table |5] The final state with three mini-scars of the form 
5-7-5 has the lowest energy of all final states and provides a measure of E*^nneaied- 

It is also informative to track the position dependence of the interstitial energy 
after it relaxes. Numerical measurements of local relaxed energy for the interstitial 
as a function of initial location are presented in Table |6j In most cases the initial /a 
complexion undergoes transitions to more stable interstitial configurations, except 



the configuration that starts from the very center location C in Fig. 19 For the 
initial conditions, C, C", 0, and 3, I3 is the most stable. Note that the interstitial 
created at the exact center C has the lowest defect energy, while the one nearest 
to the disclination from the location 10 requires the largest energy for interstitial 
defect formation. 




Two-Dimensional Matter: Order, Curvature and Defects 



45 




Figure 23. (Colo r on line) The rotational motion of an interstitial configuration (created with the initial 
location in Fig. |19|l mediated by the transition: /s — > 74 — > /s. 



Table 6. The energy of interstitial defects created at different 
initial positions within the spherical crystal. The energies shown 
are for a power law potential with 7 — 6. 
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Figure 24. (a) Vector order on the sphere with its 2 vortex defects; (b) Nematic order on the sphere with 
its 4 type 1/2 disclination defects. Taken from I128| . 



3.6. Spherical Nematics 

The case of nematic order on the sphere is also of great interest from both a funda- 
mental and materials science viewpoint |60| I128j . The free energy for nematogens 
on the sphere is given by 
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where j3ij is the geodesic angle between defects i and j with winding numbers nj and 
rij, and E{R) is the defect self-energy. Nematic spheres might be made by coating 
a droplet with gemini lipids, ABA triblock copolymers or nanorods. Microfluidic 
techniques for creating a thin spherical shell of liquid crystal in a double emulsion 
have been explored in [ 129] I130j . Although the tetrahedral array of four type-1/2 
nematic disclinations might be unstable to Boojum formation (type-1 defects) for 
spheres immersed in a nematic bulk fluid [131J, Huber and Stark |132] have shown 
that nematic wetting for spheres immersed in an isotropic fluid is a promising 
alternative route to a spherical nematic. 

An important motivation for exploring nematic order on spherical topologies is 
the search for superatoms of low valency. We have seen that spherical crystals have 
12 distinguished regions, each region being marked by the presence of either an 
isolated 5-disclination or an extended defect array like a scar. By functionalizing 
these regions one could hope to engineer 12-valent supermolecules. In particular one 
might be able to take advantage of the unusual presence of 5-7 pairs (dislocations) 
in scars to convert scars to sticky zones of the spherical superatom. Twelve is very 
likely too high a coordination number, however, to produce viable molecules or bulk 
materials. How can we change the coordination number of spherical superatoms? 
The solution lies in changing the local order on the sphere. An ordered state on 
the sphere with p-fold symmetry will have 2p net topological defects by Eq. ([s]). 
The crystalline case discussed so far corresponds to p = 6. Vector-like order on the 
sphere corresponds to p = 1 and the defects correspond to index-one vortices - 
one can also visualize the source and sink of fluid flow on a sphere [see Fig. [24)^ a)] 
or imagine the base and crown bald spots of hair combed on a sphere. Thus a 
2-sphere decorated with magnetic (or any other vector) field lines is a candidate 
for a valence-2 superatom. Consider now nematic order: p = 2. In this case the 
vortex configurations above each break up into 2 type-1/2 nematic disclinations, 
leading to a total of 4 defects located at the corners of a tetrahedron (in the one 
Frank constant approximation) [59l (see Fig. 24 ^b)). This raises the possibility of 
making tetravalent superatoms with sp^ type bonding 



The precise geometry of the defect structure is important for determining what 
kind of chemical linkers could be attached to the defects to facilitate the creation 
of large scale structures like a diamond-type lattice of linked tetravalent colloids. 
Diamond lattices are desirable for photonic crystals (optical analogs of semicon- 
ductors) since they are known to possess a large photonic band gap |133] . 

Although the elementary type-1/2 disclinations favor a tetrahedral geometry if 
splay and bend deformations are degenerate in energy this changes when elastic 
anisotropy is present. In the limit of pure splay or pure bend the energy of a -|-1 
defect becomes degenerate with two 1/2 defects. Thus two -|-1 defects sitting at 
the north and south pole can be separated be cutting the sphere on a great circle 
joining the two defects and rotating by an arbitrary angle a. This gives rise to a 
set of four 1/2 disclinations lying on the same great circle and therefore coplanar 
in the embedding space |134] . Thus the tetrahedral configuration crosses over to a 
great-circle configuration as one increases the anisotropy of the elastic constants. 
This explicitly demonstrates that defect positions can be controlled by varying 
the elastic anisotropy. Thus the structure of directional bonding for the associated 
supermolecules is controllable. 
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Figure 25. (Color online) (Left) the ground state configuration of 1082 nematic rods of aspect ratio 15 on 
a unit sphere. (Right) the exact location of the four disclinations lying near one great circle. From |134| . 

4. Crystalline order on surfaces with variable Gaussian curvature and 
boundary 

4.1. Introduction 

Under specific experimental conditions amphiphilic molecules in solution, such as 
lipids or amphiphilic block copolymers, self-assemble in a spectacular variety of 
shapes including spherical and cylindrical micelles, vesicles and lamellae, together 
with more complex geometries such as uni- and multilamellar vesicles, onion vesi- 
cles, toroidal and cage-shaped micelles. The exact shape and size of these structures 
has been observed to depend on both molecular details, such as molecular size, the 
hydrophilic/hydrophobic ratio and molecular stiffness, and collective parameters, 
such as the concentration or the diffusivity of the molecules. 

As noted previously, amphiphilic membranes can exist in different thermody- 
namic states depending on the degree of orientational and positional order and 
their precise molecular constituents. The L^J^ and Pp phases, which occur in lipid 
membranes featuring the phosphatidylcholine (PC) group, have been found, in 
particular, to exhibit in-plane orientational correlations extending over 200 A, one 
order of magnitude larger than the typical spacing between PC groups [251 1135] . 
Because of thermal fluctuations, such membranes are typically corrugated, possess- 
ing spatially inhomogeneous Gaussian and mean curvature. Furthermore it is likely 
for lipid membranes to have pores providing a passage from the exterior to the in- 
terior of a vesicle. It is therefore natural to ask how variable Gaussian curvature 
and the presence of one or more boundaries affects the phenomenology reviewed 
thus far for the sphere. 

In this section we discuss crystalline order for two important surfaces with vari- 
able Gaussian curvature and (possibly) boundary, namely the bumpy surface that 
is obtained by revolving the graph of a Gaussian function around its symmetry 
axis (i.e. a Gaussian hump) and the paraboloid of revolution. The former can be 
thought as a gentle deformation of a plane and thus can serve as a playground to 
analyze the onset of behavior not found in planar systems; the latter is possibly 
the simplest two-dimensional Riemannian surface having variable Gaussian curva- 
ture and boundary and provides a setting that is simple enough to carry out a full 
analytical treatment and analyze also large curvature regimes. Furthermore, since 
paraboloidal shapes naturally occur across the air /liquid interface of a fluid placed 
in a rotating cylindrical vessel, a direct physical realization of paraboloidal crystals 
can be constructed by assembling monodisperse objects on the surface of a rotating 



liquid. In Sec. 4.5 we review a simple experiment by Bowick et al [136] in which 
such a macroscopic model for a paraboloidal crystal is constructed by assembling 
a two-dimensional soap bubble "raft" on the air /liquid interface of a water-soap 
solution, thus extending the classic work of Bragg and Nye on planar bubble rafts. 
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Purely orientational order on the Gaussian bump has been recently reviewed by 
Turner et al [631 and won't be discussed here. 



4.2. Surfaces of revolution and conformal mapping 

Before analyzing the ground state structure of a crystalline paraboloid and Gaus- 
sian bump, it is useful to review the geometry of surfaces of revolution. The notion 
of conformal mapping of Riemannian surfaces will also be used in the following 
and will be briefly reviewed in this section with special attention to its application 
to the calculation of Green functions on simply connected 2— manifolds. 

A surface of revolution M is a surface obtained by revolving a two-dimensional 
curve around an axis. The resulting surface therefore always has azimuthal sym- 
metry. The standard parametrization of a surface of revolution is: 

X = ^(r) cos (p 

V = Or) sin 4> (107) 
z = rj{r) 

where r G [0, ii] (with R possibly infinite) and (j) = [0,27r). The metric of the 



surface (107) is given by: 

ds" = [{i'f + {r]'f]dr^ + fdcP^ , (108) 

where the prime indicate a partial derivative with respect to r. The Gaussian 
curvature is a function of r only and is given by: 

^^ .w-rv) (^09) 

^(^/2 + ^/2)2 • i^uyj 

In the presence of a boundary dM the condition ^ for the total topological charge 
of any triangulation on M reads: 

VdM Vm 

Q = ^(4-c,) + ^(6-q) =6x (110) 

i=l 1=1 

where Vqm and Vm are the number of vertices on the boundary and the interior 
of the manifold respectively (with V = Vqm + Vm the total number of vertices). 
Qi,dM = 4 — Cj is the topological charge of a vertex of coordination number q located 
on the boundary, where the coordination number of a perfect triangular lattice is 
four rather than six. A surface of revolution with boundary dM = {r = R} x [0, 2tt] 
is homeomorphic to a disk {g = and h = 1) and has therefore X = 1 and total 
topological charge Q = 6. This topological constraint can be satisfied, for instance, 
by placing six isolated 3— fold disclinations along the boundary and keeping the 
interior of the surface defect-free or by placing a 5— fold disclination in the interior 
and the remaining five 3— fold disclinations along the boundary. 

Let us consider now a generic Riemannian surface M and two curves on S in- 
tersecting at some point xq. The angle between the two intersecting curves is, by 
definition, the angle between the tangents to these curves at xq. A mapping of a 
portion S" of a surface onto a portion S* is called conformal (or angle-preserving) 
if the angle of intersection of every arbitrary pair of intersecting arcs on 5* is the 
same as that of the corresponding inverse images on S at the corresponding point 
(see for example |137j ). It is not difficult to prove that a mapping from a portion 
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5" of a surface onto a portion S* is conformal if and only if, when on S and S* the 
same coordinate systems have been introduced, the coefficients g*^ and gij of the 
metric tensor of S* and S are related by: 

9ij = w{x)gij , (111) 



with w a positive function of the coordinates x = {x^,x'^). Eq. ( 111 ) implies indeed 
that the angle between any pair of intersecting curves is the same in 5* and S. 
Isometrics are a special case of conformal mappings where w = 1 and the mapping 
is both distance and angle-preserving. 

A special type of conformal mapping is that of a portion 5 of a surface into a 
plane. This may be accomplished by introducing a set of coordinates u = (m^,m^) 
such that: 

ds"^ = w{u)[{du^f + [du'^f] . (112) 



Coordinates ((u}^v?) satisfying Eq. (112) are called isothermal (or conformal). In 



general, any simply connected Riemannian manifold with a C°°— smooth metric 
ds"^ can be equipped with a set of local isothermal coordinates. This important 
result can be stated by saying that any simply connected Riemannian manifold is 
locally conformally equivalent to a planar domain in two-dimensions. In conformal 
coordinates, the Gaussian curvature reads: 

^^_2A[logM.)]_ 
w{u) 

Conformal mapping is the fundamental tool behind the celebrated uniformization 
theorem according to which, every simply connected Riemannian surface is confor- 
mally equivalent to the unit disk, the complex plane or the Riemann sphere. This 
theorem, first proved by Koebe and Poincare independently in 1907, extends the 
Riemann mapping theorem for simply connected domains in the complex plane 
to all simply connected Riemannian surfaces and provides an insightful classifica- 
tion scheme. Its formidable power lies in the fact that the mapping that allows 
one to transform a generic surface into a simpler "irreducible" one is not an arbi- 
trary homeomorphism but is conformal, and thus preserves part of the geometrical 
structure of the original manifold. In the following we will see how this feature has 
important consequences in the elastic theory of defects on curved surfaces. The uni- 
formization of Riemannian surfaces is historically the first example of geometriza- 
tion. In the case of manifolds of higher Hausdorff dimension, and 3-manifolds in 
particular, the latter program has become, following Thurston, one of the most 
challenging and fascinating chapters of modern geometry. 

A bounded Gaussian bump and a paraboloid of revolution are both conformally 
equivalent to the unit disk D of the complex plane. Calling z = ge^'^, the new metric 
will be: 

ds^ = w{z){dQ^ + Q^d<t?) . (114) 



The conformal factor w can be found by equating the metrics ( 108 ) and ( 114 ). This 
yields: 



w{p) 



i{ry^ 



(115) 
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with Q and r related by the differential equation 



dg 
dr 



, 



(116) 



whose solution is given by: 



exp< ± dr 



(117) 



The sign of the exponent and the integration constant in Eq. (117) can be tuned 
to obtain the desired scale and direction of the conformal map. 

The calculation of the Green function of the Laplace and biharmonic operator 
is considerably simplified once a surface of revolution has been endowed with a 
local system of isothermal coordinates. In this case it is easy to show the Laplace- 
Beltrami operator Ag takes the form: 



A. 



where A is now the Laplacian in the Euclidean metric tensor: 



lee 



1, 



7e0 = 0, 



(118) 



(119) 



with determinant 7. Since the determinant of the metric tensor is rescaled (y^ 
w y/^) under conformal mapping, the Laplace and biharmonic equation for the 
Green function become: 

AGl{z,0=5{zX) (120a) 
Aw-^AG2l{z, C) = 8{z, C) (120b) 
where 5{z^ C) is the standard delta function at the point z = Q ol the unit disk. 



Eq. (|120a|) is now the standard Laplace-Green equation. Its associated Dirichlet 

z-C 



problem has the familiar solution: 

Gl[zX) 



1 
2^ 



log 



l-zC 



121 



Eq. (120b) is known as the weighted biharmonic Green equation. The uniqueness of 



its solution requires imposing both Dirichlet and Neumann boundary conditions: 

(122) 



rG2L(z,C)=0 zE^D 
19„(,)G2l(^,C) =0ze5D 



where denotes the derivative with respect to the variable z along the normal 
direction at 9D. Its solution can be expressed in integral form as 



G2l{z,0 = j d\GL{z,a)[GL{aX)-H{a,0] , 



(123) 



where H(a, Q is a harmonic kernel that enforces the Neumann condition. Such a 
function depends on the form of the conformal weight w. For radial weights w = 
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i/;o(|zp), such as those obtained by conformany mapping a surface of revolution, 
i?(c7, C) has been calculated explicitly by Shimorin |138] : 

H{aX) = -2 fjj^ dswo{s)k(^^C^) (124) 



where: 



n>0 ^" n<0 "^l*^! 



and the coefficients Cn are given by: 



Cn = 2 [ dtewoit) . (126) 
Jo 



4.3. Crystalline order on the Gaussian bump 

The most natural way of introducing a non-vanishing Gaussian curvature on an 
initially flat medium is to gently deform the medium at one point in such a way 
that the curvature introduced by this deformation dies off at infinity. If the initial 
planar domain is the entire Euclidean plane R^, the Euler characteristic is zero and, 
independently of the value of the Gaussian curvature, an ordered phase embedded 
on it is not topologically required to contain defects. Such a construction is clearly 
ideal to detect the onset of structural behavior not occurring in the ground state of a 
planar system such as the appearance of dislocations and disclinations. Vitelli et al 
|139| 1140] analyzed crystalline and p— atic order on the bumpy surface obtained by 
revolving the graph of a Gaussian function about its symmetry axis. The resulting 
Gaussian hump has parametrization: 

' X = r cos ^ 

< y = rsiii(j) (127) 
z = /iexp(-^) 

with r G [0, R] and (j) G [0, 2tt). In this parametrization the metric tensor gij (with 
determinant g) and the Gaussian curvature are given by: 

grr = i{r), 5r</. = 0, g^^ = r"^ ^ 

rim V rl)' 
where a = h/vQ \s the aspect ratio of the bump and l{r) is given by 

2 2 2 

^(r) = l + ^e"%. (129) 

^0 



(128a) 
(128b) 
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It is instructive to verify that the Euler characteristic x vanishes when the boundary 
radius R is set to infinity. Employing the Gauss-Bonnet theorem one has 



X 



[ dr^K+^ I dsKg , (130) 
Jo 27r Jc^ 



where Kg, the geodesic curvature of the circular boundary Cr of radius i?, is given 

by 

^ (131) 



If the bump is unbounded the second term in Eq. ( 130 ) disappears and, in the limit 
R ^ oo, one has 

POO 

Jo 

which implies x = 0- For finite values of R, on the other hand, the first integral in 



Eq. ( |130D gives: 



[ dry/g 
Jo 



K = 1 



so that the terms proportional to 1/ y^i{R) cancel each other, yielding x = 1- 
Because the infinite bump has x = disclinations must appear in pairs and dislo- 
cations must have total Burgers vector zero: 

i i 

Vitelli et al showed that topological defects appear in the ground state when the 
aspect ratio a exceeds a critical value Oc. Because of the topological constraint, de- 
fects appear initially in the form of a pair of unbound dislocations, roughly located 
in the region where K = Upon increasing the aspect ratio, more dislocations 
appear. This mechanism clearly resembles the defect proliferation that occurs in 
two-dimensional melting and suggests an interpretation of the curvature as a local 
effective temperature. 

At the onset of defect proliferation, inter-defect interactions are negligible com- 
pared to the interactions with the "smeared out" topological charge associated with 
the curvature of the underlying medium. The dislocation unbinding mechanism is 
then described by a non-local function of the Gaussian curvature representing the 



defect-curvature interaction part of the elastic energy of Eq. ( 54 ) : 



Pint = Y j d^xr,{x) j d^yG2L{x,y)K{y) =Y j d^xrj{x)ip{x) , (132) 

where Fint is the curvature-defect interaction part of the elastic energy and ip{x) 
can be interpreted as a geometric potential associated with the Gaussian curvature 
of the embedding manifold: 



A^ip{x) = K{x) 



(133) 
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Figure 26. (Color online) Curvature-defect interaction energy of an isolated disclination for a Gaussian 
bump of aspect ratio a = 0.5. Open symbols represent the data from a numerical mi nimi zation of a 
fixed connectivity harmonic model. The lower and upper branch are obtained from Eq. ( |136| by setting 
9 = ±tt/2 and letting A equal 4 (blue curve) and 8 (red curve). [Courtesy of V. Vitelli, University of 
Pennsylvania, Philadelphia, PA]. 



Thus 



^{x) = J (fy Gl{x, y) [Gl^V, z)K{z) + U{y)] , (134) 

where U{y) is a harmonic function which enforces the boundary conditions. The 
Laplacian Green function has the form (121) with the conformal distance p = |z| 
given here by 



= ^ exp 



(135) 



For a single dislocation of Burgers vector b, Fmt has the form [139J : 



int 



1 e 
-hba^Y sind) 




(136) 



where A = r/ro and A = R/tq. The first term in Eq. (136) corresponds to the 
R ^ oo geometric potential, while the second is a finite size correction arising from 



a circular boundary of radius R. A plot of the function ( 136 ) is shown in Fig. 26 The 
profile of the function Fint can be undestood by regarding a dislocation as bound 
pair of disclinations of opposite topological charge. Each disclination interacts with 
a potential of the form (133). For small r, positive (negative) disclinations are 



attracted (repelled) by the center of the bump. As a consequence disclinations 
experience a force which increases linearly with r. Thus if the positive disclination 
in the dipole is closer to the top, it will experience a force that is opposite and 
slightly less than that acting on the negative disclination that is further from the 
top. As a result an effective "tidal" force will push the dislocation downhill. For 
large r, however, the geometric potential saturates and the attractive force exerted 
on the positive disclination takes over and drags the dislocation toward the center 
of the bump. The minimum of the elastic energy corresponding to the equilibrium 
between these two competing forces is obtained for A « 1.1. The origin of these 
forces is the Peach-Koehler force fk = ^kj^iCTij acting on a dislocation of Burgers 
vector hi in an external stress field aij. In the case of a two-dimensional crystal 
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Figure 27. (Color o nline ) Dislocation unbinding, (a) Critical aspect ratio Qc as a function of ro/fe. The 
theoretical estimate (|137|l is plotted versus ro/fe for core energies Ej, = (dashed line) and Be = ClS^y 
(solid line). Circles are obtained from a numerical minimization of a fixed connectivity harmonic model. 
On right the logarithm of the numerically calculated strain energy for a (c) defect-free configuration on a 
bump of ro = lOf) and a = 0.7 > ac and (d) the defective configuration shown in the inset. [Courtesy of 
V. Vitelli, University of Pennsylvania, Philadelphia, PA]. 



on a substrate with preexisting Gaussian curvature, dislocations couple with the 
internal stress due to the curvature of the medium, with similar effects. 

Unbound dislocations occur in the ground state of a crystalline Gaussian bump 
when the aspect ratio a exceeds a critical value. For nearly flat landscapes the 
energetic cost of a pair of unbound dislocations is larger than that arising from 
the distortion of the medium and the system is favored to be defect free. Upon 
increasing the aspect ratio, however, the resulting elastic strain can be partially 
relieved by introducing a pair of dislocations with equal and opposite Burgers 
vector. The transition occurs when the energy gain from placing each dislocation 
in the minimum of the potential energy Fint outweights the total work needed to 
tear them apart plus the core energies 2E(i. The resulting critical aspect ratio is 
given by: 



h , /ro 



(137) 



where b' = {b/2)e ^'^E^_/{Yb'^) _ rpj^g number of unbound dislocations increases with 
increasing aspect ratio and elementary 5—7 dislocations cluster in more complicated 



structures with zero net Burgers vector. Fig. 27 shows a plot of the critical aspect 



ratio Oc as the function of the dimensionless parameter ro/6 from Ref. |139j . as 
well as density plots of the strain energy corresponding to a defect-free bump (c) 
and a representative configuration featuring two unbound dislocations (b). 

The interaction between defects and curvature also has some remarkable sta- 
bilizing effects on defect dynamics. Dislocation dynamics consists of two distinct 
processes: glide and climb. Glide is motion along the direction of the Burgers vector. 
It requires only local rearrangement of atoms and thus has a very low activation 
energy, making it dominant at low temperatures. Climb consists of motion in the 
direction perpendicular to the Burgers vector. It requires diffusion of vacancies 
and interstitials and is usually suppressed therefore relative to glide. As a conse- 
quence of the underlying curvature, however, a gliding dislocation experiences a 
"recalling" force f recall ~ kd\y\i where y is the transverse displacement and kd is 
a position-dependent effective spring constant. At leading order in y and a, the 
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Figure 28. (Color online) Dislocation glide, (a) Recalling potential ^k^iy'^ acting on dislocations gliding 
in the direction y (expressed in units of the lattice spacing a) for ro/a = 10, A = 0.5 (at y = 0), A = 8 and 
a = 0.1 (dark blue), 0.3 (blue), 0.5 (orange) and 0.7 (red). The energy is scaled by 10~*. Effective spring 
constant for tq = 10a and a = 0.5. The ordinate axis is scaled by 10~^. Filled and empty circles represent 
numerical data obtained from a fixed connectivity harmonic model. [Courtesy of V. Vitelli, University of 
Pennsylvania, Philadelphia, PA]. 



latter is given by: 



kd{r) 



ba^Y 
4ro 



1-(1 + A2)e- 
A3 



(138) 



This effective recalling force is not due to any external field nor to the interaction 
of dislocations with other defects, but exclusively to the coupling of the gliding 
dislocation with the curvature of the substrate. As expected, the effective spring 
constant (138) vanishes for planar crystals (i.e. a — > 0). Since Yb'^ can be hundreds 
of ksT, at finite temperature the harmonic potential ^kdy"^ associated with the 
recalling force raises the activation energy of thermally induced dislocation glide. 



The harmonic recalling potential and the effective spring constant (138) are plot- 



ted in Fig. 28 together with numerical data obtained from a fixed connectivity 
harmonic model [139j. 

Spatial curvature also provides an effective potential in the thermal diffusion of 
interstitials and vacancies. These can be constructed by grouping three dislocation 
dipoles. The elastic energy associated with a single interstitial/vacancy can be 
derived from Eq. (132) in the form 



1, 



Fintix) « -YnV{x 



(139) 



where V = Aip and 0, is the area excess or deficit associated with the defects. 
As a result, interstitials tend to climb to the top of the bump while vacancies are 
pushed into the flat regions. Like a disclination of positive topological charge, an 
interstitial is attracted to regions of positive Gaussian curvature while a vacancy 
is attracted to regions of negative Gaussian curvature. 

The crystalline structure arising on a Gaussian bump has been investigated nu- 
merically by Hexemer et al |141] . By means of smart Monte Carlo (SMC) sim- 
ulations, the authors analyzed the ground state configuration of a system of 
point-like particles (with N up to 780) interacting on a Gaussian bump of variable 
aspect ratio with a Yukawa potential of the form U{r) = exp(— Kr)/r, with r the 
Euclidean distance between two particles in R^. This potential approximately de- 
scribes the interaction between charged particles in solution with counter-ions, with 
being the Deybe-Hiickel screening length. Starting from an initially defect-free 
lattice and relaxing it with order 10^ SMC iterations, Hexemer et al observed the 



appearance of defects for increasing aspect ratio. Fig. 29 shows the arrangement of 
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Figure 29. (Color online) Delaunay triangulation of 780 particles on a Gaussian bump, (a) and (b) Top 
view of the lattice containing two dislocations of opposite Burgers vector on a bump of aspect ratio 
a = 0.82. (c) and (d) Three dislocations arranged around a bump of aspect ratio a = 0.95. (e) and (f) A 
more complex arrangements of dislocations on a bump of aspect ratio a = 1.58. [Courtesy of Alexander 
Hexemer, UC Santa Barbara, Santa Barbara, CA. Currently at Lawrence Berkeley National Lab., Berkeley, 
CA]. 



780 particles for three different aspect ratios. As predicted by the elastic theory, 
the proliferation of defects starts with the appearance of a pair of isolated dislo- 
cations (Figs. 29 i and b). They are sitting at A = 1.46 from the center and are 



rotated by 180°, giving rise to a configuration with total Burgers vector close to 



zero. Upon increasing the aspect ratio, a third dislocation is observed (Fig. 29 



and d). The three dislocations are rotated by 120° with respect to each other. The 



green and blue circles in Fig. 29 ; and d mark the region of zero Gaussian curvature 



and the minimum of the geometric potential at 1.1 tq. As shown in the figure, for 
a = 0.95 two of the three dislocation are not sitting at 1.1 tq as one might expect. 
One is deep inside the area of positive Gaussian curvature while the other is in the 
region of negative Gaussian curvature. Hexemer et al attributed this discrepancy 
to finite-size effects. 

and f show the lowest energy state on a bump of aspect ratio a 



Fig. 29 



1.58. 

In this situation dislocations appear arranged in the form of grain boundaries. This 
phenomenon has a simple interpretation. The total elastic energy of a collection of 
defects on a curved surface consists of three parts: the curvature-defect interaction 
discussed above, which tends to localize the defects where the geometric poten- 
tial is minimal, a defect-defect interaction which is repulsive for like-sign defects 
and attractive for defects of opposite sign and an additional contribution due to 
the curvature alone. For a small number of dislocations, inter-defect interactions 
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Figure 30. Total number of disclinations in the ground state as a function of the aspect ratio a. [Courtesy 
of Alexander Hexemer, UC Santa Barbara, Santa Barbara, CA (currently at Lawrence Berkeley National 
Lab., Berkeley, CA).] 



are negligible compared with curvature-defect interactions and dislocations tend 
to gather where Fint is minimal. As shown from the data in Fig. [30] however, the 
total number of dislocations in the bump is linearly proportional to the aspect 
ratio. When the total number of dislocations on the bump exceeds some threshold 
value the repulsive interaction takes over and dislocations start spreading. Dislo- 
cations with parallel (antiparallel) Burgers vector are organized in such a way as 
to maximize (minimize) their reciprocal distance while still taking advantage of 
the Gaussian curvature screening. The competition leads to the appearance of the 



y— shaped grain boundaries shown in Fig. 29 The branching allows the disloca- 
tions to lower their potential energy by keeping a large number of dislocations close 
to the 1.1 ro circle. 




4.4. Paraboloidal crystals 

A paraboloid of revolution is a simple surface with both variable Gaussian curvature 
and boundary. Its standard parametrization is given by 



(140) 



with r £ [0,R] and (p = [0,27r]. Here h is the height of the paraboloid and R the 
maximum radius. In the following we will call k = 2h/B? the normal curvature of 
the paraboloid at the origin. The metric tensor gij and the Gaussian curvature are 
given respectively by: 

g„ = 1 + K^r'^, gr(f> = 0, g<t,<i, = r^, (141a) 

The problem of finding the optimal arrangement of disclinations in the ground 
state of paraboloidal crystals has been considered by the authors of this review 
article |142| 1143] . As in the case of a Gaussian bump, when the maximal Gaussian 
curvature exceeds some critical value (depending on k and R) defects proliferate 
in an initially defect-free configuration. Unlike the Gaussian bump, however, the 
Gaussian curvature on a paraboloid is strictly positive and only vanishes at infinity. 
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Positive isolated disclinations are thus energetically preferred to dislocation dipoles 
and the defect proliferation mechanism consists of a "migration" of one of the 
six topologically required +l-disclinations from the boundary to the origin of the 
paraboloid. Dislocations, on the other hand, appear in the system clustered in the 
form of grain boundary scars in the high density regime as in spherical crystals. 
Since the Gaussian curvature is not constant throughout the manifold, however, 
this transition is preceded by a regime in which isolated disclinations and scars 
coexist in the crystal. 

Let p{x) be the effective topological charge density of a system of N disclinations 
on a background of Gaussian curvature K{x): 



N 



(142) 



i=l 



If free boundary conditions are choosen, the elastic energy (55) can be written as 

Fei = ^J d^xr\x) , (143) 

where T(x) = Axix), and x(^) satisfies the inhomogeneous biharmonic equation 

A\{x) = Yp{x) , (144) 

with boundary conditions 



X{x) = X £ dM 
u'Vixix) X G dM 



(145) 



where is the ith component of the tangent vector perpendicular to the bound- 
ary. In the parametrization (140), the normal vector v is simply given by gr/loA, 



with Qr the basis vector associated with the coordinate r. The stress function r(a;) 
can thus be expressed as 



EN 

Y 



(fyGL{x,y)p{y) + U{x). 



(146) 



Gl{x, y) is given in (121 ) and U{x) is a harmonic function on the paraboloid that 



enforces the Neumann boundary conditions in Eq. ( 145 ) . The conformal distance 
on the paraboloid is given by 



Q{r) = A 



re 



1 + Vl + 



(147) 



Using Eqns. (147) and (121) in Eq. (143), the elastic energy of a collection of 



N disclinations in a paraboloidal crystal can be expressed as: 
^ = ? E "i^Gdx, X,) - T,{\x\) + U{x) , 



(148) 



i=l 



where the first term represents the bare contribution of the defects to the energy 
density and the second corresponds to the screening effect of the Gaussian curva- 
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Figure 31. (Color online) The function ri,(r) and Ui^^u for different values of k. 

ture. Explicitly: 



rs(|a^|)=log 



1 + Vl + K^r"^ 



(149) 



where r = Iccl and 



a 



1 + + K?R'^ 
exp fVl + K2i?2 



(150) 



is a normalization constant depending on boundary radius R and the ratio k. Fig. 
313, shows a plot of the screening function Ts{r) for different values of k G [1,2]. 
As expected, the contribution due to Gaussian curvature is maximal at the origin 
of the paraboloid and drops to zero at the boundary. 

The calculation of the harmonic function U{x) requires a little more effort. If the 
crystal was defect-free (or populated by a perfectly isotropic distribution of defects) 
the function U {x) would be azimuthally symmetric and constant on the boundary. 
By the maximum principle of harmonic functions, U{x) would then be constant 
on the whole manifold and depend only on k and the radius R: U{x) = U^^r- This 
constant can be determined by integrating Ax{x) = T{x) and imposing the second 



boundary condition in Eq. (145). This gives: 



27r 
~A 



dr^/gTsir), 



(151) 



where A is the area of the paraboloid: 



A 



2tt 
3^ 



'1 + k'R') 



(152) 



As shown in Figure [31p , the value of C/k,_r quickly approaches the linear regime as 
the size of the radius increases: 



K,R 



1 1 

-- nR H 

4 3 



(153) 



Then, for a defect-free configuration, the contribution of the boundary to the en- 
ergy density is a constant offset that persists even for large radii. In the presence 
of disclinations, on the other hand, the function x(^) is no longer expected to 
be azimuthally symmetric and the harmonic function U{x) will not be constant 
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throughout the paraboloid. In this case U{x) can be expressed in the integral form 

U{x) = - I d^yH{x,y)p{y), (154) 



where H{x,y) is the harmonic kernel (124). In the regime of large core energies 



Fc ^ Fgi , the creation of defects is strongly penalized and the lattice necessarily has 
the minimum number of disclinations allowed by the topology of the paraboloidal 
substrate. From symmetry considerations, we might expect the optimal distribution 
of defects to consist of b +1— disclinations arranged along the boundary at the base 
vertices of a 6-gonal pyramid and a 6-fold apex (of topological charge qq = 6 — b) at 
the origin. The homogeneous boundary conditions adopted require the first term 



in Eq. ( 148 ) to vanish at the boundary. In the minimal energy configuration then, 



the system has the freedom to tune the total number of defects along the boundary 



to minimize the elastic energy Eq. ( 143 ) for any given value of the ratio k. This 



behavior is exclusive to manifolds with boundary and doesn't have any counterpart 
in crystals on compact surfaces like the sphere and the torus. 

We will label a pyramidal configuration by Yj,, where b denotes the number of 
base +1— disclinations. The coordinates {r,(p) of the vertices are given by 

Yf. |(0,any), (155) 

Using the Euler theorem one can show that it is possible to construct infinite 
families of polyhedra with the symmetry group C^y from the pyramidal backbone 
Yf,. The number of vertices is given by: 

V = lbn{n + l) + 1, (156) 

where n is a positive integer which represents the number of edges (not necessarily 
of the same length) of the polyhedron which separates two neighboring disclina- 
tions. In the following we will refer to these polyhedra with the symbol Yfj^^- Fig- 
32 illustrates two l^^n lattices for the cases 6 = 4 and n = 7 (with V = 113), and 
6 = 5 and n = 10 {V = 276). By a numerical minimization of the energy Eq. ( 143| ) 



one can establish that the Y^ are indeed equilibrium configurations for 6 G [3,5], 
for some range of the parameters k and R. The cases of 6 = 5, 6 are particularly 
significant because they are characterized by an equal number of defects {N = 6) 
of the same topological charge {q = 1). The two configurations will be associated 
therefore with the same core energy Fc and this introduces the possibility of a 
structural transition between Y^ and Ig governed by the curvature ratio k and the 
boundary radius R. For fixed R and small values of k, the 6— fold symmetric con- 



figuration Yq is the global minimum of the free energy Eq. (143). For k larger than 
some critical value Kc{R), however, the Yq crystal becomes unstable with respect to 
the 5— fold symmetric configuration Y^. A numerical calculation of the intersection 
point between the elastic energies of Y^ and Yq for different values of n and R allow 



us to construct the phase diagram shown in Fig. 33, The word "phase" in this 
context refers to the symmetry of the ground state configuration as a function of 
the geometrical system parameters k and R. In principle, if we keep increasing the 
curvature we might expect the crystal to undergo a further transition to the I4 
phase. In this case, however, the core energy will also increase by a factor 4/3 and 
so this is not generally possible in the regime in which Fc ^ Fd. For intermediate 
regimes (i.e. Fc ~ F^i), Y^ — > I4 and I4 I3 transitions are also possible. The 
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Figure 32. (Color online) Two examples of Y{,^„ triangulations of the paraboloid (5^,7 on the top and 
^5,10 on the bottom). Plaquettes with disclinations are highlighted in red, for +1— disclinations, and green 
for +2— disclinations. 




0.5 



1.5 



R 



Figure 33. (Color online) Phase diagram in the large core energy regime. For small k the lattice preserves 
the 6— fold rotational symmetry of the flat case. As the curvature at the origin increases the system 
undergoes a transition to the Y5 phase. 



critical value of the parameters k and i?, however, is not universal and will depend 
on the precise values of the core energy and the Young modulus. 



When the core energy is small, the elastic energy Eq. ( 143 ) can be lowered by 



creating additional defects. Let us assume that a fivefold disclination is sitting at 
the point xq = (ro, (/)o)- We can introduce a notion of distance on the paraboloid by 
setting up a system of geodesic polar coordinates (s, (^) with origin at xq. We expect 
that the stress introduced by the defect is controlled by an effective disclination 
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charge inside a circular domain of geodesic radius L: 

<leff = q- dip I ds^/gK{s,ip), (157) 
Jo Jo 

where q = vr/S is the charge of the isolated defect and the integral measures the 
screening due to the total Gaussian curvature within the domain. The metric tensor 
and the Gaussian curvature of a generic Riemannian manifold can be expressed in 
geodesic polar coordinates in the form (see for example Do Carmo |il44'j ): 

gss = 1, gsip = 0, g^^ = G, (158a) 

«2 77^ 

K{s,v) = -^, (158b) 



where G = g<p ■ gip- Furthermore, an expansion of the metric around the origin 
(0, ip) yields: 

VG = s-IKos^ + o{s^). 



For small distance from the origin, Eq. (157) becomes: 

f27r 



q^ff = q+ dip dsdl^/G (159) 
Jo Jo 

= q-TTKoL'^ + o{L^). (160) 



The right hand side of Eq. ( 160 ) is a very general expression for the effective discli- 
nation charge at small distance and doesn't depend on the embedding manifold. If 
a grain boundary is radiating from the original disclination, we expect the spacing 
between consecutive dislocations to scale like a/qeff, with a the lattice spacing 
[75] . When qeff — > O"*" the dislocation spacing diverges and the grain boundary 
terminates. Since the Gaussian curvature is not constant, the choice of the origin 
(i.e. the position of the central disclination along the grain boundary) affects the 
evaluation of qeff - One can identify upper and lower bounds by observing that: 

maxi^(r) = K{0) = k^, (161a) 

r 

minK(r) = i^(i?) = (^^^J^- (161b) 

Unlike the case of surfaces of constant Gaussian curvature, the phase diagram for 
paraboloidal crystals consists of three regions separated by the curves: 

KqL^ = 1 Ko = i^^in, Kmax. (162) 

When L — L{K^ia) O"*", the effective disclination charge goes to zero and the 
distance between two consecutive dislocations diverges at any point. On the other 
hand if L — L(i^jnax) 0~, the disclination charge will prefer to be delocalized in 
the form of grain boundary scars. For L(i('max) < L < L(A'jnin) the paraboloid will 
be equipped with regions where the Gaussian curvature is high enough to support 
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Figure 34. (Color online) Defect phase diagram for a paraboloidal crystal of radius R = 1. The two 
phase boundaries that separate the isolated disclinations (ID) regime from th e co existence regime and 
the coexistence regime from the scar phase correspond to the solutions of Eq. | |163| for Kq = K^^in and 
Ko = ifmax, respectively. 



the existence of isolated disclinations as well as regions where the screening due to 
the curvature is no longer sufficient and the proliferation of grain boundary scars 
is energetically favored. This leads to a three region phase diagram in which the 
regime of isolated disclinations is separated from the delocalized regime of scars 
by a novel phase in which both isolated disclinations and scars coexist in different 
parts of the paraboloid according to the magnitude of the Gaussian curvature. 
It is useful to measure the distance Lc in terms of the lattice spacing a and 



rephrase Eq. (162) as a condition on a (or equivalently on the number of vertices 
V). To do this we note that in order for the domain Cl to completely screen the 
topological charge of the shortest scar possible (i.e. 5 — 7 — 5), the geodesic radius L 
has to be large enough to enclose the entire length of the scar. Galling i the geodesic 
distance associated with a single lattice spacing a, we will then approximate L ~ 3i. 
This leads to the following expression for the lattice spacing a at the onset of scar 
formation: 

"'"iC-^^iTi)^ ''''' 



The lattice spacing a can be approximately expressed as a function of the number 
of vertices of the crystal by dividing the area A of the paraboloid by the area of a 
hexagonal Voronoi cell of radius a/2 with ~ V. The phase diagram arising 

from the solution of Eq. (163) is sketched in Fig. p4]for the case R = 1. The two 
phase boundaries that separate the isolated defects (ID in the plot) regime from 
the coexistence regime and this one from the grain boundaries phase correspond 
to the solutions for Kq = Kmm and Kq = i^max respectively. 

Fig. 35 shows the Voronoi lattice and the Delaunay triangulations obtained from 
a numerical minimization of a system of V classical particles (up to V = 200) 
constrained to lie on the surface of a paraboloid and interacting with a Coulomb 
potential of the form U{r) = 1/r, with r the Euclidean distance between two 
particles in R^. Simulations are carried out using a parallel implementation of Dif- 
ferential Evolution (DE) |142| 1145] . In all the systems observed disclinations always 
appear clustered in either grain boundary scars or dislocations with the exception 
of isolated +1— disclinations which appear in the bulk as expected from the curva- 
ture screening argument discussed above. The complex aggregation of defects along 
the boundary together with the presence of negatively charged clusters indicates 
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Figure 35. (Color online) Voronoi lattice and Delaunay triangulations for ten selected systems from nu- 
merical simulations with R = 1. The first row corresponds toV = 200 and k = 1.6, while the second row is 
for V = 200 and k = 0.8. In the bottom four rows V = 150 , 100 , 80 , 60 , 50 , 30 , 20 , 16 and k = 1.6. From 
[1431. 



that the effect of the boundary, in the case of relatively small systems like the ones 
simulated, is more drastic than predicted by the homogeneous boundary conditions 
Eq. (145). Even in the computationally expensive case V = 100, the distance 
between the origin and the boundary of the paraboloid is only four lattice spacings. 
In this situation we expect the distribution of particles along the boundary to play 
a major role in driving the order in the bulk. 

For larger systems, such as V = 200 (top of Fig. 35), the behavior of the parti- 
cles in the bulk is less affected by the boundary and the crystalline order reflects 
more closely the free-boundary problem previously presented. A comparison of the 



35 


in 
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K = 0.8 and V = 200, the defects are 
all localized along the boundary with the exception of one length-3 scar in the bulk 
at distance r ~ 0.63 from the center. For k = 1.6, the pattern of defects in the bulk 
is characterized by the coexistence of an isolated -|-1— disclination at the origin and 
a length-5 VF— shaped scar displaced along a parallel one lattice spacing away from 
the central disclination. Apart from the evident difficulty in comparing the struc- 
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tures of small systems with those predicted from continuum elasticity theory, this 
behavior is consistent with the simple picture sketched in the phase-diagram of Fig. 



34, The local 5— fold symmetry at the origin of the k = 1.6 configuration, compared 



with 6— fold symmetry for k = 0.8, suggests, as in the case of spherical crystals 
[M] . that the complicated structure of defect clusters appearing in large systems 
is the result of the instability of the simpler 1^^^ configurations from which they 
partially inherit their overall symmetry. A more accurate numerical verification of 
our theory remains a challenge for the future. 

The symmetry of the configurations presented in Fig. 35 deserves special atten- 
tion. As for any surface of revolution, the circular paraboloid possesses the sym- 
metry group 0(2) of all rotations about a fixed point and reflections in any axis 
through that fixed point. Any given triangulation of the paraboloid may destroy 
the full rotational symmetry completely or just partially, leaving the system in one 
of the following two subgroups: the pyramidal group Cnv or the reflection symme- 
try group Cs- In general we found the latter symmetry group for system sizes up 
ioV = 200 particles. The symmetry for larger system sizes is under investigation. 



4.5. Experimental realization of paraboloidal crystals 

Some sixty years ago Bragg and Nye used bubble rafts to model metallic crys- 
talline structures |146]. A carefully made assemblage of bubbles, floating on the 
surface of a soap solution and held together by capillary forces, forms an excel- 
lent two-dimensional replica of a crystalline solid, in which the regular triangular 
arrangement of bubbles is analogous to the close packed structure of atoms in a 
metal. Feynman considered this technique to be important enough that the famous 
Feynman lectures in physics include a reproduction of the original Bragg-Nye paper 
in its entirety jl47 ]. Bubble rafts can be made easily and inexpensively, equilibrate 
quickly, exhibit topological defects such as disclinations, dislocations and grain 
boundaries, and provide vivid images of the structure of defects. Bubble raft mod- 
els have been used to study two-dimensional polycrystalline and amorphous arrays 
|148) . nanoindentation of an initially defect-free crystal |149] . and the dynamic be- 
havior of crystals under shear |150] . In this section we review the experimental 
realization of a bubble-raft model for a paraboloidal crystal done by Bowick et al 
1 136] by assembling a single layer of millimeter-sized soap bubbles on the surface 
of a rotating liquid, thus extending the classic work of Bragg and Nye on planar 
soap bubble rafts. 

As we mentioned in the introduction, paraboloidal shapes naturally occur across 
the air/liquid interface of a fluid placed in a rotating cylindrical vessel. It is simple 
to verify that the height z of such an interface above the xy— plane of a system of 
Cartesian coordinates is given by 



^2 

z=—r^ + C, (164) 
25 



where oj is the angular velocity and g the gravitational acceleration and C = 
D — uPB?' /Ag is an integration constant following from the requirement that the 
volume of the liquid during rotation be equal to the volume of the liquid at rest {D 
is its height) |151] . Thus the normal curvature at the origin is given by k = uP' /g. 



Fig. 36 shows a computer reconstruction from Ref. |136J of two bubble rafts 
with K = 0.15 cm~^ and k = 0.32 cm^ respectively. Soap bubbles have diameter 
a = 0.84(1) mm with monodispersity Aa/a ~ 0.003. The images shown in Fig. 



36 are taken with a CCD digital camera rotating along with a cylindrical vessel 
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Figure 36. (Color online) Lateral and top view of a computer reconstruction of two paraboloidal rafts 
with Ki 0.15 cm~^ (a, b) and ft2 ~ 0.32 cm"'^ (c, d). The number of bubbles is A'^i = 3813 and 
N2 = 3299 respectively. The color scheme highlights the 5— fold (red) and 7— fold (blue) disclinations over 
6— fold coordinated bubbles (yellow). From Ref. [136) . 



of radius R = 5 cm and Delaunay triangulated to determine the lattice topology. 
The degree of order of the crystalline raft can be characterized by the translational 
and orientational correlation functions g{r) and ge^r) [llOj . The former gives the 
probability of finding a particle at distance r from a fixed particle at the origin. The 
function is normalized with the density of an equivalent homogeneous system in 
order to ensure g{r) = 1 for a system with no structure. Interactions between par- 
ticles build up correlations in their position and g{r) exhibits decaying oscillations, 
asymptotically approaching one. For a two-dimensional solid with a triangular lat- 
tice structure the radial correlation function is expected to exhibit sharp peaks 
at r/a = y/n'^ + nm + m? = 1, \/3j 2, 2\/3 • • •, while the amplitude of the peaks 
decays algebraically as r"*^ with r/ = 1/3 (dashed line in Fig. 37). Within the pre- 
cision of the data, the positional order of the paraboloidal crystals assembled with 
the bubble raft model reflects this behavior, although with more accurate measure- 
ments one might see dependence of the exponent rj on the curvature (because of 
the proliferation of defects with curvature). 

The orientational correlation function gQ{r) is calculated as the average of the 
product {ijj{Qi)'ijj*{r)) of the hexatic order parameter over the whole sample. For each 
bubble (labeled j) that has two or more neighbors, '4>j{'r') = {^/Zj) "^/.L^ exjp{6i9jk), 
where Zj is the number of neighbors of i and Oj^ is the angle between the j — k 
bonds and a reference axis. One expects gei^r) to decay exponentially in a disordered 
phase, algebraically in a hexatic phase and to approach a finite asymptote in the 
case of a crystalline solid. In the systems studied gei^r) approaches 0.8 within 5~6 
lattice spacings. 

Of particular interest is the structure of the grain boundaries appearing in the 
paraboloidal lattice for different values of the curvature parameter k. Grain bound- 
aries form in the bubble array during the growing process as a consequence of 
geometrical frustration. As noted, any triangular lattice confined in a simply con- 
nected region with the topology of the disk is required to have a net disclination 
charge Q = 6. In the absence of curvature, however, the elastic stress due to an 
isolated disclination is extremely high and defects are energetically favored to clus- 
ter in the form of a grain boundary consisting of one- dimensional arrays of tightly 
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Figure 37. (Color online) Translational and orientational correlation functions (g and t/gi respectively) for 
rafts with (a,b) k ^ 0.32 cm"!, a = (0.8410± 0.0025) mm, and (c,d) k 0.15 cm"!, a = (0.9071 ± 0.0037) 
mm. All the curves are plotted as functions of r/a, where r is the planar distance from the center and a 
is the bubble radius. The envelope for the crystalline solid decays algebraically (dashed line), while the 
orientational correlation function approaches the constant value 0.8. From Ref. I136| . 




5 mm 

Figure 38. (Color online) An enlarged view of the terminating grain boundary scar shown in Fig |36| 1 
for a system with large Gaussian curvature. The scar starts from the circular perimeter of the vessel 
and terminates roughly in the center carrying a net +1 topological charge. The image of the bubbles (left) 
shows that they may deform slightly to better fill space, whereas the computer reconstruction of the lattice 
(right) uses perfect spheres of uniform size. From Ref. 11361 . 

bound (5, 7)— fold disclination pairs. In a planar confined system, grain boundaries 
typically span the entire length of the crystal. In the curved case, however, they can 
appear in the form of scars carrying a net +1 topological charge and terminating 
in the bulk of the crystal. 
Prominent examples of grain boundaries are visible in the two lattices shown in 

0.15 cm-^l 



Fig. 



36 



For a gently curved paraboloid (with k ~ 0.15 cm^ ), grain boundaries 
form long (possibly branched) chains running from one side to the other and passing 
through the center. As the curvature of the paraboloid is increased, however, this 



long grain boundary is observed to terminate in the center (see Fig. 361; a close- 



up version of this image is seen in Fig. 38). For R = 5 cm, the elastic theory of 



defects predicts a structural transition at Kc = 0.27 cm~ in the limit of large 
core energies. In this limit the creation of defects is strongly penalized and the 
lattice has the minimum number of disclinations required by the topology of the 
embedding surface. In a low curvature paraboloid (k < Kc) these disclinations 
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are preferentially located along the boundary to reduce the elastic energy of the 
system. When the aspect ratio of the paraboloid exceeds a critical value Kc{R), 
however, the curvature at the origin is enough to support the existence of a 5— fold 
disclination and the system undergoes a structural transition. In the limit of large 
core energies, when only six disclinations are present, such a transition implies a 
change from the Cq^ to the C^^ rotational symmetry group. 

Together with the theoretical argument reviewed in the previous section, these 
experimental observations point to the following mechanism for scar nucleation in 
a paraboloidal crystal. In the regime in which the creation of defects is energeti- 
cally inexpensive, geometrical frustration due to the confinement of the lattice in a 
simply connected region is responsible for the formation of a long side-to-side grain 
boundary. When the curvature of the paraboloid exceeds a critical value, depen- 
dent on the radius of the circular boundary, the existence of a -|-1 disclination near 
the center is energetically favored. Such a disclination serves as a nucleation site for 
5-7 dislocations and the side-to-side grain boundary is replaced by a terminating 
center-to-side scar. 



5. Crystalline and p— atic order in toroidal geometries 
5.1. Introduction 

Circa twenty years afer their first observation in partially polymerized diacetylenic 
phospholipid membranes by Mutz and Bensimon |152j . self-assembled toroidal ag- 
gregates are now considered the progenitors of a magnificent cornucopia of complex 
structures, also featuring branched network and micellar surfaces of high genus. The 
existence of such complex structures has become an experimental fact thanks to 
the enormous work in the past decade on the study of self-assembly of amphiphilic 
compounds such as lipids, surfactants and amphiphilic block copolymers. After 
the work of Eisenberg and coworkers |153| 11541 1155] it became clear that block 
copolymers in particular afford access to a variety of complex structures spanning 
an unexpectedly vast range of topologies and geometries. More recently, several 
experimental studies have been performed on di- and triblock copolymers with 
the intent of unraveling the origin of morphological complexity in copolymer sur- 
factants and some possible pathways for micelle and vesicle formation have been 
proposed. 

Jain and Bates, for instance, reported the formation of several non-simply- 
connected micellar structures from the self-assembly of diblock copolymer poly(l,2- 



butadiene-6-ethylene oxide) (PB-PEO) [156" (see Fig. 39). These polymers self- 
assemble in y— shaped junctions and these form the building blocks for more com- 
plex micellar structures via re-assembly. Later, Pochan et al. jl57j found that al- 
most all of the microstructures assembled from poly(acrylic acid-6-methyl acrylate- 
6-styrene) (PAA99-PMA73-PS66) triblock copolymers are ringlike or toroidal mi- 
celles. The route to toroidal micelles in block copolymers, however, is believed to 
be different from that to more complex network structures suggested by Jain and 
Bates, since residual y— shaped aggregates are very rarely found in the sample. 
According to the authors of Ref. [ISTj, the formation of toroidal micelles has to 
be attributed primarly to the collapse of cylindrical micelles. On the other hand, 
mere end-to-end connection of cylindrical micelles doesn't appear to be the exclu- 
sive ring-forming mechanism in triblock systems since the average circumference 
of the self-assembled tori appears smaller than the contour length of the average 
cylindrical micelles. It seems to be a more complicated process, possibly involving 
interactions and exchange of matter between neighboring cylinders. 
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Figure 39. Complex non-simply connected structures from self-assembly of diblock copolymer poly(l,2- 
butadiene-b-ethylene oxide) (PB-PEO). From Ref. [156]. 



Toroidal micelles have also been observed in recent experiments by Kim et al 
|158j from the self-assembly of amphiphilic dumbbell molecules based on a aro- 
matic rod segment that is grafted by hydrophilic polyether dendrons at one end 
and hydrophobic branches at the other end. Molecular dumbbells dissolved in a 
selective solvent self-assemble in an aggregate structure due to their amphiphilic 
character. This process has been observed to yield coexisting spherical and open- 
ended cylindrical micelles. These structures, however, change slowly over the course 
of a week to toroidal micelles which thus appear more stable. The formation of ring- 
like structures was explained in this case as result of the coalescence of spherical 
micelles occurring to reduce the contact between hydrophobic segments and water 
molecules. 

An alternative pathway for the self-assembly of toroidal structures in copolymer 
solutions has been proposed and numerically tested by He and Schimd |159] . In this 
pathway, the micelles do not coalesce, but simply grow by attracting copolymers 
from the solution. Once a critical micelle size is exceeded, copolymers start to flip- 
flop in such a way that the micelle core becomes itself solvent-philic (semi-vesicle 
state). Finally the solvent diffuses inside the core and the semi- vesicle swells into 
a vesicle. 

The existence of toroidal aggregates in amphiphilic compounds was suggested 
theoretically for smectic-C (SmC) membranes, based on the argument that ori- 
entational order is not frustrated on a surface with zero Euler characteristic and 
thus a toroidal topology may be energetically favoured in p—ai\c membranes with 
large order parameter coupling compared with the bending rigidity |59| 1160] . This 
hypothesis was tested by Evans in 1995 with the result that a toroidal topology 
is indeed prefered over the spherical one for wide range of geometrical and me- 
chanical parameters in defect-free p— atic tori. The role of disclinations in p— atic 
tori was investigated systematically by Bowick, Nelson and Travesset nine years 
after the original work of Evans with the conclusion that, even if not required by 
topological constraints, unbound disclinations can be energetically convenient even 
in very dense systems ^161] . The precise number of unbound disclinations is con- 
trolled primarly by the aspect ratio of the torus. The existence of defects in the 
ground state of an ordered phase on the torus is crucial for toroidal crystals where 
it leads to some unique structural features as well as a spectacular example of a 
curvature-driven transition to a disordered, liquid-like, state in the limit the aspect 
ratio approaches to one |162| 1163] . In the remainder of this section we will review 
these three examples of order on embedded tori. 
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Figure 40. (Color online) The standard parametrization of a circular torus of radii and i?,2. 



5.2. Geometry of the torus 

The standard two-dimensional axisymmetric torus embedded in is obtained by 
revolving a circle of radius R2 about a coplanar axis located at distance Ri > R2 
from its center. Choosing the symmetry axis as the z direction of a Cartesian frame, 
a convenient parametrization can be found in terms of the polar angle ip along the 
revolving circle and the azimuthal angle (p on the xy— plane: 

X = + R2 cos ip) cos (p 

y = {Ri + R2 cos ip) sin (j) , (165) 
z = R2 sin Ip 

where ip ,(p & [— vr, vr). These coordinates satisfy the Cartesian equation 

(^Ri - ^Jx"^ + + z'^ = Rl (166) 

In this parametrization the metric tensor gij and Gaussian curvature K are given 
by 

g^^ = RI, 9^4, = ^, 94,4, = {Ri + R2 cos V')^ , (167) 

cos tp 

^ " R2iRi + R2cosij) ■ 

The Gaussian curvature is therefore positive on the outside of the torus, negative 
on the inside and zero along the two circles of radius Ri at ip = ib7r/2 (see Fig. 



40). Moreover K is maximally positive along the external equator at ip = and 
maximally negative on the internal equator at ^p = itTr. The Gauss-Bonnet theorem 
requires the total topological charge and the integrated Gaussian curvature to be 



zero on the torus: 

TV 

^qk= / (fxK{x) = . (169) 
k=i 

A global measure of the curvature of the embedded torus is provided by the aspect 
ratio r = Ri/R2- Our discussion will be limited to the case r > 1. A "fat" torus 
with r = 1 is obtained by taking Ri = R2 and is characterized by a singularity 
in the Gaussian curvature along the internal equator at ■0 = ibvr. The case r < 
corresponds to a self-intersecting torus whose symmetry axis lies in the interior of 
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the revolving circle. The "skinny" torus limit, r — > oo, can be obtained either by 
taking i?2 ^ with finite Ri or by letting R2 stay finite and taking i?i — > 00. The 
Gaussian curvature diverges in the first case and goes to zero in the second. Neither 
case, however, reflects a real physical situation since the area A = {2tt)'^ R1R2 of 
the torus becomes zero or infinite respectively. 
Non-axisymmetric tori can be generated in various ways by deforming the surface 



of revolution (165). A special and physically relevant transformation consists of an 
inversion wrt the unit sphere, a translation about a vector (3 and a second inversion. 
Thus a point R on the surface is mapped onto: 

This map is conformal and transforms an axisymmetric shape to a non-axisymmet- 
ric one if the vector /3 is not parallel to the symmetry axis |164] . As /3 is increased 
in magnitude, the asymmetry increases until, in the limit, the torus becomes a 
perfect sphere with an infinitesimal handle. The physical importance of mapping 
(170) relies on the fact that it leaves the bending energy k J cPxH^ invariant with 



significant consequences for the equilibrium shape and stability of fluid toroidal 
membranes. 

As in any closed manifold the traditional Green-Laplace equation doesn't have 
a solution on the torus. An alternative Green function, can be obtained from the 
modified equation: 

AgGL{x,y) = Sg{x,y)-A-^ (171) 

As usual the calculation of the Green function can be simplified considerably by 
conformally mapping the torus to a domain of the Euclidean plane via a suitable 
system of isothermal coordinates. Intuitively the torus is conformally equivalent 
to a rectangular domain described by a system of Cartesian coordinates. To make 
this explicit, one can equate the metric of the torus in the coordinates {ip, (/>) to a 
conformally Euclidean metric in the coordinates ^): 

ds'^ = Rldip'^ + {Ri + R2 cos -ipfdct? = w {df + drf) , 

where w is a positive conformal factor. Taking r] = (j) and w = -|- -R2 cos V')^, 
the coordinate ^ is determined by the differential equation: 



d^ 1 



dtp r + cos ip ' 



(172) 



where r = R1/R2, the aspect ratio of the torus, may be taken greater than or equal 
to one without loss of generality. Choosing the plus sign and integrating both sides 



of Eq. (172) we find: 



Jq r + cos ip' 



Taking ip G [— vr, vr), the integral (173) yields: 



(. = K arctan lu tan 
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where 



2 /r - 1 
K= , t^=W • (174) 

In the transformed coordinate system r]) the modified Green-Laplace equation 
reads: 

AGL{x,y) = d{x,y)-j, (175) 

where A and S are now the Euchdean Laplacian and delta function. The function 
GL{x,y) can be expressed in the form: 

Gl{x, y) = Go{x, y) - {G^{x, ■ )) - (Go(- , y)) + (Go(- , • )) , 

where Go{x,y) is the Laplacian Green function on a periodic rectangle and the 
angular brackets stand for the normalized integral of the function Go{x,y) with 
respect to the dotted variable: 

{Go{x,-)) = J ^Go{x,y). (176) 
Analogously the function (Go(- , • )) is given by 

(Go(-,-)) = / ^^Go{x,y) 
and ensures the neutrality property: 

J (fx Gl{x, y) = j (fy Gl{x, y) = 0. (177) 

The modified Laplacian Green function on a periodic rectangle of edges pi and p2 
can be conveniently calculated in the form: 

o„(.,y) = ^'-MEM, (178) 

where ux is the eigenfunction of the Laplace operator with periodic boundary 
conditions: 

Auxix) = Xuxix) , (179) 

such that: 

jux{0,r)) = ux{pi,v) 
\uxi^,0) = uxi^,P2) ' 

In Cartesian coordinates the eigenfunctions are simple plane waves of the form: 
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where A„ and Hm are given by: 



27rn 27rm 
An = — lJ'm = — — n, m = 0, ±1, ±2 



Pi P2 
and the eigenvalue A is given by: 



A - -^2 _ „2 

Calling for simplicity x = {x,y) and y = (C,??), the function Gq is given by: 



:i8ii 



Go{x,y) 



E 



^i\„{x-£,)^ifj„„{y-ri) 



Summing this series eventually leads to |163] : 



Goix,y) 



log 2 1 



R o \y -V\ + T^^oi 
DVT 2 27r 



where '!9i(n|r) is the Jacobi theta function |165| 1166] defined as: 

00 

'&i{u\t) = 2g4 sin w n ~ '^^^ + ^'^") ~ ^ 

n=l 

with g = exp(i7rr), z = x + iy and C = C + 



2n\ 



(182) 



(183) 



(184) 



5.3. Defect- free p—atic textures and genera transition 

There has been much effort in recent years to shed light on possible pathways in 
the self-assembly of toroidal micelles and vesicles from block copolymer solutions 
and other amphiphilic compounds. In particular, the spontaneous formation of 
structures of genus g > 1 from preexisting spherical objects, a process appearing 
in several different scenarios proposed in the literature, has attracted the most 
attention and debate because of its exotic character. The simplest question one 
can ask in this context is whether a vesicle can be energetically favored to change 
its topology from spherical to toroidal once the total surface area, and therefore the 
number of constituent molecules in the vesicle, is specified. Evans addressed this 
problem in 1995 |160j and showed how such a transition between genera is indeed 
possible in fluid membranes and is even enhanced if the membrane is endowed 
with in-plane orientational order. Although oversimplified if compared to the great 
complexity of the real self-assembly mechanism, Evans' calculation provides insight 
into the delicate problem of stability of toroidal vesicles as well as a good starting 
point for our discussion of order on the torus. 

Let 9 be the local orientation of a p—atic director field on a surface as defined in 
Sec. 2. A standard orientational order parameter is given by scalar field: 



(185) 



where (•) deontes a thermal average. The total elastic energy of the vesicle consists 
of a pure bending term Ff,, describing the elasticity of the membrane in the liquid 
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state, and a term Fp associated with the internal p— atic order. Thus F = Fp + F^, 
with: 



Fp = J <fx (r|V^|2 + + C|(V - ipil)^^- 



(186a) 
(186b) 



where k and the bending and Gaussian rigidity respectiveljj^ and the op- 

erator V — ipQ, is obtained from the covariant derivative of a p— atic tensor order 
parameter expressed in a local orthonormal frame as described in Sec. |2.4[ with Q, 



the covariant vector of Eq. (53 ) that parametrizes the spin-connection. The free en- 



ergy ( 186a) is invariant under rotations of the local reference frame by an arbitrary 
angle x- 



^ - diX ■ 



(187a) 
(187b) 



which is a typical gauge transformation. The gradient term in Eq. (|186a|) is the 
same as in Eq. 



52 



upon identifying Ka = p^C. As observed by Park et al. |168j . 



Eqns. (186a) and (186b) closely resemble the Landau-Ginzburg Hamiltonian of a 



superconductor in an external magnetic field: 



n 



V 



2ie 

he 



H |V X A 

8^' 



(188) 



When exposed to an external magnetic field, a super-conducting material can un- 
dergo a second order mean-field transition from a metal to a super-conductor char- 
acterized by an Abrikosov lattice of vortices whose density is determined by the 
temperature and the applied magnetic field H. The magnetic field, in particular, 
is conjugate to the vortex number since: 



d^xV X A = (t>oLNy , 



(189) 



where L is the length of the sample in the direction of H and (j)o = hc/2e is the flux 
quantum. In this context the vector potential associated with the applied magnetic 
field is replaced by the covariant vector Q, whose curl is the Gaussian curvature. 
Morover, on a closed surface: 



J d^xV xfl = J 



d^xK 



27r 
P 



N 



i=l 



Thus a atic phase on a closed 2— manifold is analogous to an Abrikosov phase 
with a fixed total vorticity rather than a fixed applied magnetic field. 



^The expression of the bending energy used in Eq. ( |186b| l is that originally gave by Helfrich for a membrane 
with zero pre-existing curvature 11671 . Often the equivalent expression J d?x {^kH^ + KgK) is found in the 
literature. In this case, however, the mean curvature H is defined as the sum of the principal curvatures 
rather than their average. 
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The transition to an ordered phase in controlled by the parameter r in Eq. ( 186a 



For T above a critical value Tc {jc = on a flat surface) ,■0 = and the system is in 
the isotropic phase. In this regime the only contribution to the energy is determined 
by the shape of the vesicle as expressed by the bending term (186b). For r < Tc, on 



the other hand, the gauge symmetry is spontaneously broken and ■i/' 7^ 0, with the 
exception of a number of isolated defective points. Within a meanfield approach, 
thermal fluctuations can be neglected and the preferred conflguration of the system 
corresponds to the minimum of the free energy. A ground state conflguration for 
the atic order parameter can be found by rewriting the gradient energy term in 
Eq. (186a) as 



with Dk = dk — ip^k, and expressing ip in the basis of eigenfunctions of the operator 
in parentheses: 



an^n 



(190) 



with ifn satisfying the Hermitean equation 

1 



-D [yjg Dk) = 

9 ^ 



with the normalization condition 



d^x <flfr. 



A 

4tt 



(191) 



(192) 



The complex coefficients are then determined by minimization of the free en- 
ergy F. The eigenfunctions are sometimes referred to as "Landau levels" and are 
typically degenerate. If ■0 is expanded in this complete set of eigenfunctions then, 
in meanfleld, the partition function is dominated by conflgurations involving only 
Landau levels with the lowest eigenvalue Aq. Taking lowest levels only, the atic 



free energy (186a) becomes 



d^x 



(193) 



To flnd the lowest energy conflguration of the fleld ip one now has to solve the 



eigenvalue problem (191) to determine the lowest eigenvalue Aq and minimize the 
free energy with respect to the parameters of the linear combination ( |190 ) . Notice 
that the mean fleld transition has been lowered from Tc = to Tc = —^t^XqC/A. 
Although Ao = on a flat surface, it is flnite and positive-deflnite on a surface with 
flnite Gaussian curvature [160'. 

For a sphere of unit radius parametrized in standard spherical coordinates {6, (p), 
the eigenvalue problem ( 191 ) is solved hy Xq = p and ipn of the form |169| 
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for integer values of n between —p and p. These functions have 2p zeros, corre- 
sponding to topological defects of unit charge and winding number 1/p. Because 



of the lowest Landau level approximation used to derive Eq. (193), eigenfunctions 



(194) are degenerate to quadratic order in ijj. This implies the freedom to place 
defects anywhere on the sphere with no additional energy cost of order ilP' . The ■0^ 
term, on the other hand, lifts this degeneracy and make the defects repel. 



For axisymmetric tori the eigenvalue problem (191) was solved numerically by 



Evans with the result shown in Fig. 41 Eigenvalues are plotted as a function 
of the aspect ratio r of the torus for various p— atic textures labelled pn with 
< n < p. Eigenfunctions have 4n zeros corresponding to 2n disclinations of 
topological charge g = 1 (winding number 1 /n) and 2n disclinations of topological 
charge q = —1 (winding number — 1/n). Except for vector order {p = 1), defective 
configurations are energetically favored for small aspect ratios and the number of 
disclination pairs increases with r as a consequence of the larger (in magnitude) 
Gaussian curvature. 



Carrying out the integrals in Eq. (186a) and (186b) in the lowest Landau level 



approximation, the total elastic free energy of spherical and toroidal vesicles can 
be written as 



i^sphere = 47r(2K + Kg) - 27r — 



27rV 



2^2 



27r 



C72 
uA 



uA J i/sphcrc(p) 
JtornsiP,r) 



(195a) 
(195b) 



where J is the minimal value of the integral of |V'mf|^ and 



Yll=-p (^nVn sphere 



a^(fp + a-tp-n torus 
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Figure 42. (Color online) Phase diagram for vesicles of constant area of genus zero and one with intrinsic 
vector (top left), nematic (top right), tctradic (bottom left) and hexatic (bottom right) order. The position 
of the transition line depends on the value of Kg/n. Toroidal vesicles that exist in the shaded region are 
non-axisymmetric. Data are taken from |160| . 



is the linear combination of eigenfuctions that are degenerate to quadratic order 
in ^l^. The second term in Eqs. (195a) and (195b) disappears in the isotropic phase 
when r > —^^ttXqC/A due to the intrinsic orientational order on the manifold. 
Even in this case, Eqs. ( |195a ) and (195b) reveal the existence of a transition line 
between spherical and toroidal shape. In the absence of in-plane orientational order 
the bending energy term in Eq. ( 195b ) is minimized by the so called Clifford torus 
with r = ^/2. In this case -Ftorus = 47r'^K, which is smaller than i^sphere = 47r(2K-|-«;g) 
when Kg/ K, > vr — 2. Toroidal vesicles with r = \/2 are therefore energetically fa- 
vored in the fluid phase for large values of Kg/n. This argument cannot be invoked 
to explain the complex self-assembly of toroidal structures from homogeneous so- 
lutions mentioned in the introduction, but does provide a simple (but non-trivial) 
example of a toroidal shape being energetically preferred to a spherical one. 

The critical value of Kg/n of the genera transition is lowered for vesicles possessing 
in-plane p— atic order. This is displayed in the phase diagrams of Fig. |42] for the 
case of vector {p = 1), nematic {p = 2)^ tetradic {p = 4) and hexatic {p = 6) order. 
The lines, whose position depends on Kg/n, separate spheres (above the line) from 
tori (below the line). Even for simple vector order, stable toroidal vesicles exist as 
well for Kg/n < tt — 2. For Kg/n > vr — 2, furthermore, transition lines are closed, 
with spheres on the inside and tori on the outside. In the shaded regions non- 
axisymmetric tori are favored over symmetric ones. The hexatic phase diagram 
(bottom right corner of Fig. 42 ) deserves special attention. The shaded region is 
split in two parts. That on the right contains spheres if Kg/n < tt — 2 and non- 
axisymmetric tori otherwise, while that on the left contains spheres above the 
transition line and non-axisymmetric tori below. The white region of the phase 
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diagram is also divided in two. The top left region of the diagram behaves as 
described above, while the other part is characterized by closed lines separating 
spheres (outside) and small tori (inside). These last class of tori exhibit ten pairs of 
g = =bl disclinations with positive disclinations distributed on the external equator 
of the torus and negative disclinations along the internal equator. This last feature 
is an important property sheared by toroidal objects with in-plane order and will 
be clarified in the following sections. The results reviewed here are valid within the 
mean-field approach and the lowest Landau level approximation, which both hold 
deep in the ordered phase that we are most interested in. At higher temperatures 
Evans proved the approximations are still valid away from the transition line and 
for C/A 3> n/ksT ^ 1 [160]. The latter condition is generally fulfilled by several 
systems (i.e. K/ksT = 1 — 10 for lipid bilayers). 



5.4. Defective ground states in hexatics 

Evans' analysis, summarized in the previous section, indicates that disclinations, 
even if not required by topological constraints, can nonetheless appear in the 
ground state of p— atic tori as a consequence of the coupling between in-plane 
orientational order and spatial curvature. The occurrence of defects in the ground 
state of toroidal hexatic vesicles was systematically investigated by Bowick, Nelson 
and Travesset based on the formalism outlined in Sec. 2 |161] . Before embarking 
on a detailed analysis of the elasticity of defects in hexatic tori, it is instructive to 
obtain a rough estimate of how many pairs of ibl-disclinations would be required 
to achieve perfect screening of the background topological charge associated with 
the Gaussian curvature of the torus. Consider a wedge of angular width A0 on the 
outside wall of positive Gaussian curvature. The net curvature charge associated 
with this region is 

Seff= / d<p di;^K = 2A4>. (196) 

Upon equating Se// to 2tt/6, the charge of a single disclination, one finds that A(f> = 
27r/12, independently of iii and i?2- Thus, 27r/A0 = 12 positive disclinations would 
be required to completely compensate the negative curvature of the inner wall. 
This simple argument neglects core energies and interactions between disclinations, 
effects which will cause the preferred number of defect pairs to be less than twelve. 

The total elastic energy of a toroidal hexatic vesicle containing N disclinations 
of topological charge gj (i = 1 . . . A^) takes the form 



„ 2 AT 

i<j 1=1 



27r^ K 27r^r^ ec 2 



where the first two terms represent the contributions due to the pair interaction 
between defects and the interaction of defects with the topological charge of the 



substrate associated with the Gaussian curvature. The third term in Eq. (197) 
is the spin-wave part of the frustrated hexatic energy while the last two terms 
represents the bending energy and the defect core energy respectively. The defect- 
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Figure 43. (Color online) The total energy (in units of Ka/'2 of hexatic tori of aspect ratio r = \/2 
(left) and r = 2.6926) for varying number of defects. The bending energy at fixed r is subtracted off. The 
disclination core energy is set to 0.1X4, which is 0.2 in the above units. 

defect interaction potential has been calculated in Ref. [161J: 



Q{Xi,Xj) = - — log 



27r 



i Oi—aj 



2tt yfW^ 



4^2 



T] 



(198) 



"2 _ 1 V 27r 

where a is related to the polar angle on a cross-section of the torus by: 

r cos a — 1 , , 

cosV = (199 

r — cos a 

and r] is the Dedekind eta function: 



77(r) 



e 24 



2TTinT\ 



(200) 



n=l 



The one-body interaction between defects and curvature, on the other hand, is 
expressed via the potential energy 



C{xi) = log 



1 



r — cos a,- 



(201) 



For opposite sign disclinations, the defect-defect interaction, determined by the 
function Q{xi, xj), is attractive for all separations at constant (p. If only this term 
were present, the attraction would bring both charges as close as possible, binding 
all disclinations into dipoles which have a higher energy than a defect-free con- 
figuration. Thus if no other terms were present, the ground state configuration 
would be defect free. The defect-curvature interaction term C{xi), however, favors 
the appearance of additional defects. This term acts like an electric field pulling 
the positive (negative) disclinations into regions of positive (negative) Gaussian 
curvature. 



The elastic energy (197) was analyzed numerically in Ref. [16V for various aspect 
ratios and defect core energies. A plot of the energy of configurations obtained by 
placing a ring of {N/2) (-l-l)-disclinations equally spaced along the same parallel on 
the outside of the torus and a second ring of {N/2) (— l)-disclinations on the inside. 
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is shown in Fig. 43 for various N and variable angular separation A-i/^ between the 
rings. The disclination core energy is taken to be Cc = cKa, with c a numerical 
constant taken equal to 0.1 in the plots. The energy at the maximum separation 
(A'!/; = tt) first decreases and then increases with A^. The optimal number of defect 
pairs N/2 is 6 and 7 (the latter curve is not shown in the plot) for r = \/2 and 
r = 2.6926 respectively. 

The existence of defects in the ground state of a toroidal vesicle with hexatic order 
is also affected by the total number of molecules N forming the hexatic phase. In 
Ref. |161j it was found that defects disappear in the limit ^ cxd, but are present 
for numbers of molecules as large as = 10^". This number is several orders of 
magnitude larger than the typical number of molecules of biological vesicles such 
as red-blood cells (i.e. N ~ 10^). To make this estimate, one considers a pair 
of opposite sign disclinations that have been pulled apart from the circle of zero 
Gaussian curvature to the equators in the regions of like sign Gaussian curvature; 
thus V+ = and V- = t^- Upon approximating the defect-pair energy by its flat 
space value, the total energy reads: 



TT 



i^^log 



r + l 
r — 1 



+ j^Ka log — 
18 V OQ 



(202) 



where oq is the lattice spacing. Eq. (202) changes sign when 
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ao 
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Taking 
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leads to the conclusion that defects are favored for: 
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(204) 



where the last identity has been obtained by taking Ec/Ka = 0.1. This result 
establishes that defects are present in the ground state of a hexatic tours for any 
fixed number of molecules provided the torus is sufficiently fat. For the energetically 
favored Clifford torus with r = ^/2, Eq. ( 204 ) predicts a critical number of molecules 
to be order Af ~ 10^°. 

In absence of defects the ratio between the hexatic stiffness Ka and the bending 
rigidity k dictates the optimal shape of the toroidal vesicle. Taking qi = in Eq. 



(197) one obtains in this case: 
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2tt^Ka 
+ Vr'^-1 



+ K- 



27r^r 



2^2 



(205) 



The first term represents the energetic cost associated with the distortion of the 
hexatic director field due to the Gaussian curvature alone. In the limit of large and 
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which provides a compelling example of the interplay between order and geometry: 
if the stiffness associated with the intrinsic hexatic order is much smaller than 
the bending rigidity, the Clifford torus is the optimal geometry; if on the other 
hand, the hexatic stiffness dominates, then a thin torus, similar to a bicycle tire, 
is optimal. 



5.5. Toroidal crystals 

Crystalline assemblages of identical sub-units packed together and elastically bent 
in the form of a torus have been found in the past ten years in a variety of systems 
of surprisingly different nature, such as viral capsids, self- assembled monolayers 
and carbon nanomaterials. In the introduction we mentioned the self-assembly of 
toroidal micelles and vesicles from homogeneous solutions of amphiphilic molecules 
such as oligomers of aromatic compounds or block copolymers. Toroidal geometries 
also occur in microbiology in the viral capsid of the coronavirus torovirus |170] . 
The torovirus is an RNA viral package of maximal diameter between 120 and 140 
nm and is surrounded, as other coronaviridae, by a double wreath/ring of cladding 
proteins. 

Carbon nanotori form another fascinating and technologically promising class 
of toroidal crystals [49j with remarkable magnetic and electronic properties. The 
interplay between the ballistic motion of the tt electrons and the geometry of the 
embedding torus leads to a rich variety of quantum mechanical properties includ- 
ing Pauli paramagnetism |171] and Aharonov-Bohm oscillations in the magnetiza- 
tion |172j . Ring closure of carbon nanotubes by chemical methods |173] suggest 
that nanotubes may be more flexible than at first thought and provides another 
technique of constructing carbon tori. In this section we review some recent de- 
velopments in the study of the geometry and the elasticity of toroidal crystals. 
Additional details can be found in Refs. |162| 1163] . 

5.5.1. Geometry of toroidal polyhedra 

Toroidal crystals with local 6— fold order can be described as toroidal polyhedra 
with all triangular faces. After the discovery of carbon nanotubes in 1991 and the 
subsequent theoretical construction (later followed by the experimental observa- 
tion) of graphitic tori, this class of polyhedra has drawn considerable attention 
and many possible tessellations of the circular torus have been proposed by the 
community [Ei [ESI ETSl II771 [Ml [ESI EBOl EED EEa [I83l EEl] . Here we review 
the construction of a defect-free triangulated torus and we show how the most 
symmetric defective triangulations can be generally grouped into two fundamental 
classes corresponding to symmetry groups and D^d respectively. 

The structure of a triangulated cylinder can be specified by a pair of triangular 
lattice vectors c and t, called the chiral and translation vector respectively, which 
together define how the planar lattice is rolled up. In the canonical basis ai = (1,0) 

and a2 = (2, ^), the vector c has the form 



c = nai -\- ma2 n, m £ Z 



(206) 
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Figure 44. (Color online) Construction of a defect-free triangulation of the torus. On top planar map of 
the triangulated torus corresponding to the (n, m, Z) configuration (6,3,1). On the bottom (6,3,6) chiral 
torus. The edges of each one of the six tubular segments has been highlighted in red. 

The translation vector t, on the other hand, can be expressed as an integer multiple 



t = let I e 



(207) 



of the shortest lattice vector et perpendicular to c. The vector is readily found 
to be of the form 



(n + 2m)ai — (2n + m)a2 
(n + 2m : 2n + m) 



where (a : b) denotes the greatest common divisor of a and b and enforces the 
minimal length. The three-dimensional structure of the torus is obtained by con- 
necting the edge OT of the rectangle in Fig. ^ (left) to (TC and OC to (TT. The 



edge OT is then mapped to the external equator of the torus while the edge OC 
to the (p = meridian. The resultant toroidal lattice has characteristic chirality 
related to the initial choice of the vector c. In the nanotubes literature armchair 
refers to the lattice obtained by choosing n = m, zigzag to that obtained for m = 
and chiral to all other lattices. An example of an (n, m, I) chiral torus is shown in 
Fig. 44 (right) for the case n = 6, m = 3 and I = 6. The chirality is extremely 
important in graphitic carbon nanotube or nanotori, where it determines whether 
the electronic behavior of the system is metallic or semiconducting. 

By Euler's theorem one can prove that the number of triangular faces F and the 
number of vertices V a triangular toroidal lattice is given by 



V 



2^ ■ 



Denoting Aji the area of the rectangle with edges c and t and At the area of a 
fundamental equilateral triangle, the number of vertices of a defect-free toroidal 
triangulation is then: 



V 



Aji 21 (n^ -I- nm + m^) 
2At ~ {n + 2m:2n + m) ' 



(208) 



The planar construction reviewed above allows only lattices with an even number 
of vertices. Defect-free toroidal deltahedra with an odd number of vertices are 
also possible and their construction is generally achieved by assembling congruent 
octahedral building blocks. We refer the reader to Ref. |185] for an comprehensive 
review of the topic. 
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The embedding of an equal number of pentagonal and heptagonal disclinations 
in the hexagonal network was first proposed by Dunlap in 1992 as a possible way 
to incorporate positive and negative Gaussian curvature into the cylindrical geom- 
etry of carbon tubules |174] . According to the Dunlap construction the necessary 
curvature is incorporated by the insertion of "knees" (straight cylindrical sections 
of the same diameter joined with a kink) in correspondence with each pentagon- 
heptagon pair arising from the junction of tubular segments of different chirality 
(see Fig. 46). In particular, a junction between a (ra, 0) and a (m, m) tube can be 



Figure 45. (Color online) Voronoi lattices of a TPn prismatic (top) and TAn antiprismatic (bottom) 
toroids with R\ = 1 and R2 = 0.3. 



obtained by placing a 7— fold disclination along the internal equator of the torus 
and a 5— fold disclination along the external equator. Since the radii of the two 
segments of a junction are different by construction, the values of n and m are 
commonly chosen to minimize the ratio |c(„^o)l/|c(m,m)l = n/y/2im. By repeating 
the 5 — 7 construction periodically it is possible to construct an infinite number of 
toroidal lattices with an even number of disclinations pairs and dihedral symmetry 



group Dnh (where 2n is the total number of 5 — 7 pairs. Fig. 45 top and Fig. 47). The 
structure of the lattice is described by the alternation of two motifs with crystalline 
axes mutually rotated by 30° as a consequence of the connecting disclination. One 
of the fundamental aspects of Dunlap's construction is that all the disclinations are 
aligned along the two equators of the torus where the like-sign Gaussian curvature 
is maximal. As we will see below, this feature makes these arrangements optimal 
in releasing the elastic stress due to curvature. Other defective triangulations with 
Dnh symmetry group have been proposed in Ref. |163] and won't be discussed here. 

Another class of crystalline tori with dihedral antiprismatic symmetry Dnd was 
initially proposed by Itoh et al |177| I178| 1179] shortly after Dunlap. Aimed at 




Figure 46. (Color online) Dunlap knees obtained by joining two straight tubular segments with {n, 0) and 
(m,m) chirality. [Courtesy of A. A. Lucas and A. Fonseca, Facultes Universitaries Notre-Dame de la Paix, 
Namur, Belgium!. 
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Figure 47. (Color online) Five-fold polygonal torus obtained by joining (5,5) and (9,0) tubular segments 
via ten pairs of 5 — 7 rings. This structure was originally proposed by the authors of Ref. 11861 11871 as a 
possible low-strain configuration for carbon nanotori. [Courtesy of A. A. Lucas and A. Fonseca, Facultes 
Universitaries Notre-Dame de la Paix, Namur, Belgium]. 

reproducing a structure similar to the Cgo fullerene, Itoh's original construction 
implied ten disclination pairs and the point group D^^. In contrast to Dunlap tori, 
disclinations are never aligned along the equators in antiprismatic tori, instead be- 
ing staggered at some angular distance 5tjj from the equatorial plane. Hereafter we 
will use the symbols TPn and TAn to refer to toroidal lattices with 2n disclination 
pairs and symmetry group Dnh and Dnd respectively. 




(a) (b) 

Figure 48. (a) Unit cell for toroidal antiprisms. 5— fold vertices are circled and 7— fold vertices are boxed, 
(b) Unit cell of a D^^ torus in the Berger-Avron construction. The graph consists of four generation of tiles 
and the internal equator of the torus is mapped into the horizontal line passing to the mid-point between 
the 6th and the 7th vertex. 



A systematic construction of defected triangulations of the torus can be achieved 
in the context of planar graphs |188| I189| 1190] . A topological embedding of a graph 
in a two-dimensional manifold corresponds to a triangulation of the manifold if 
each region of the graph is bounded by exactly three vertices and three edges, and 
any two regions have either one common vertex or one common edge or no com- 
mon elements of the graph. The simplest example of toroidal polyhedra with Dnd 
symmetry group, featuring only 5— fold and 7— fold vertices, can be constructed by 
repeating n times the unit cell of Fig. [48^ . These toroidal antiprisms have V = 4n 
vertices and can be obtained equivalently from the edge skeleton of a n— fold an- 
tiprism by attaching at each of the base edges a pentagonal pyramid and by closing 
the upper part of the polyhedron with n additional triangles. By counting the faces 
one finds F = 5n + 2n + n = 8n from which V = 4n. The simplest polyhedron 
of this family has V = 12 and D^d symmetry group (see top left of Fig. 49) and 
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corresponds to the "drilled icosahedron" obtained by removing two parallel faces of 
an icosahedron and connecting the corresponding edges with the six lateral faces 
of an antiprism with triangular base (i.e. a prolate octahedron). Starting from 
this family of toroidal antiprisms a number of associated triangulations having the 
same defect structure can be obtained by geometrical transformations such as the 
Goldberg inclusion [3^ I19H 1192] . Such transformations, popularized by Caspar 
and Klug for the construction of the icosadeltahedral structure of spherical viruses 
|39] , consist in partitioning each triangular face of the original graph into smaller 
triangular faces in such a way that old vertices preserve their valence and new 
vertices have valence six. The partition is obtained by specifying two integer num- 
bers (-L, M) which define how the original vertices of each triangle are connected 
by the new edges so that the total number of vertices is increased by a factor 
r = L2 + LM + M2. 

A general classification scheme for Dnd symmetric tori was provided by Berger 
and Avron |189] 1190] in 1995. Their scheme is based on the construction of unit 
graphs comprising triangular tiles of different generations. In each generation, tiles 
are scaled in length by a factor l/\/2 with respect to the previous generation. This 
rescaling approximates the non-uniformity of the metric of a circular torus. 

In the past few years, alternative constructions of triangulated tori have been 
proposed as well as novel geometrical and graph-theoretical methods to express 
the coordinates of their three-dimensional structures (see for example Kirby |180) , 
Laszlo at al jl81| 1182] . Diudea et al |183] 1184] ). Here we choose to focus on the 
defect structure associated with the two most important class TPn and TAn with 
groups Dnh and Dnd- 
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5.5.2. Elasticity of defects on the torus 

The total free energy of a toroidal crystal with disclinations can be expressed 
as usual as 



1 r ^ 

i=i 



(209) 



where Y is the two dimensional Young modulus and 

N 



(210) 



fc=i 



The defect part of the stress function r(a;) can be found by integrating the Green 



function (183) and takes the form 
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where Li2 is the usual Eulerian dilogarithm and 



a 



\/ — 1 — r . 



(212) 



The function Ts{x) representing the stress field due to the Gaussian curvature of 
the torus, on the other hand, is given by 



r.(a;) 
Y 



log 



r + \/r^ — 1 
2(r + cos tp) 



+ 



r — \fr^ 
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(213) 



A derivation of the functions V^ix^Xk) and V a{x) is reported in Ref. |163j and 
won't be repeated here. 



To analyze the elastic free energy (209) we start by considering the energies of 



two opposite sign disclinations constrained to lie on the same meridian. The elastic 



free energy of this system is shown in Fig. 50 as a function of the angular separation 
between the two disclinations. The energy is minimized for the positive (5— fold) 
disclination on the external equator (maximally positive Gaussian curvature) and 
the negative (7— fold) disclination on the internal equator (maximally negative 
Gaussian curvature). The picture emerging from this simple test case suggests that 
a good ansatz for an optimal defect pattern is a certain number p of equally spaced 
+1 disclinations on the external equator matched by the same number of equally 
spaced —1 disclinations on the internal equator. We name this configuration with 
the symbol Tp, where p stands for the total number of disclination pairs. 




l<k<p 



vr, 



2TTk 
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(214) 



i<fc<p 



where the two pairs of numbers specify the {^p, (j)) coordinates of the positive and 
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Figure 50. (Color online) Elastic energy of a 5 — 7 disclination dipole constrained to lie on the same 
meridian, as a function of the angular separation. In the inset, illustration of a circular torus of radii i?i > 
R2- Regions of positive and negative Gaussian curvature have been shaded in red and blue respectively. 

negative disclinations respectively. A comparison of the energy of different Tp con- 
figurations, as a function of aspect ratio and disclination core energy, is summarized 



in the phase diagram of Fig. 51 We stress here that only Tp configurations with p 



even have an embedding on the torus corresponding to lattices of the TP| class. 
Nevertheless a comparison with p— odd configurations can provide additional in- 
formation on the stability of p— even lattices. For small core energies, moreover, 
thermally excited configurations with a large number of defects and similar p— polar 
distributions of topological charge are expected to exhibit an elastic energy com- 
parable in magnitude with that of these minimal constructions. The defect core 
energy has been expressed here in the form 



Fr 



2p 



2pec 



(215) 



The core energy Cc of a single disclination depends on the details of the crystal- 
forming material and the corresponding microscopic interactions. A simple phe- 
nomenological argument (see for example Ref. |193] ) gives 



Y 



where a the lattice spacing. Taking 
yields 
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32^ ' 



A/^V, with A the area of the torus. 
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(216) 



For a system of order F = 10^ subunits, then, the dimensionless core energy on 
the left hand side of Eq. (216) is of order 10~^. This estimate motivates our choice 



of the scale for ec/{AY) in Fig. 51 



For dimensionless core energies below 4-10~^ and aspect ratios r between 3.68 and 
10.12 the ground state structure is the TP5 lattice corresponding to a double ring 
of -|-1 and —1 disclinations distributed on the external and internal equators of the 
torus as the vertices of a regular decagon (the Tio configuration) . The TP5 lattice 
has dihedral symmetry group D^^. That this structure might represent a stable 
configuration for polygonal carbon toroids has been conjectured by the authors of 
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Figure 51. (Color online) Phase diagram for Tp configurations in the plane (r, tc/AY). For r £ [3.68, 10.12] 
and ec ~ the structure is given by a Tio configuration with symmetry group Ds^. 



Ref. [186|, based on the argument that the 36° angle arising from the insertion of 
ten pentagonal-heptagonal pairs into the lattice would optimize the geometry of a 
nanotorus consistently with the structure of the sp^ bonds of the carbon network 
(unlike the 30° angle of the 6— fold symmetric configuration originally proposed 
by Dunlap). In later molecular dynamics simulations, Han |194j, |195j found that a 
5— fold symmetric lattice, such as the one obtained from a (9,0)/(5,5) junction (see 
Fig. 47), is in fact stable for toroids with aspect ratio less then r ~ 10. The stability. 



in this case, results from the strain energy per atom being smaller than the binding 
energy of carbon atoms. Irrespective of the direct experimental observation of such 
disclinated toroidal crystals, which is still open, we show here how continuum 
elasticity predicts that a 5— fold symmetric lattice indeed constitutes a minimum 
of the elastic energy for a broad range of aspect ratios and defect core energies. 

For small aspect ratios the 5— fold symmetric configuration becomes unstable 
and is replaced by the 9— fold symmetric phase Tg. As we mentioned, however, this 
configuration doesn't correspond to a possible triangulation of the torus. It is likely 
that the ground sate in this regime consists of ten skew disclination pairs as in the 
antiprismatic TAn lattice. The latter can be described by introducing a further 
degree of freedom 5ij) representing the angular displacement of defects from the 
equatorial plane: 



TAn : 
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A comparison of the TP5 configuration and the TA5 configuration is shown in 



Fig. 52 for different values of Sif^. The intersection points of the boundary curves 



with the 5^— axis has been calculated by extrapolating the (r, 6ip) data points in 
the range 6 G [0.07, 0.8] with A{57p) = 2.57r • lO'^. For smah dip and r £ [3.3, 7.5] 
the prismatic TP5 configuration is energetically favored. For r < 3.3, however, 
the lattice undergoes a structural transition to the TA5 phase. For r > 7.5 the 
prismatic symmetry of the TP5 configuration breaks down again. In this regime, 
however, the elastic energy of both configurations rapidly rises because of the lower 
curvature and defects disappear. 

In the regime of large particle numbers, the amount of curvature required to 
screen the stress field of an isolated disclination in units of lattice spacing becomes 
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too large and disclinations are unstable to grain boundary "scars" consisting of a 
linear array of tightly bound 5 — 7 pairs radiating from an unpaired disclination 
|73| I142j . In a manifold with variable Gaussian curvature this effects leads to a 
regime of coexistence of isolated disclinations (in regions of large curvature) and 
scars. In the case of the torus the Gaussian curvature inside (IV'I > 7r/2) is always 
larger in magnitude than that outside (IV'I < 7r/2) for any aspect ratio and so we 
may expect a regime in which the negative internal curvature is still large enough to 
support the existence of isolated 7— fold disclinations, while on the exterior of the 
torus disclinations are delocalized in the form of positively charged grain boundary 
scars. 

This hypothesis can be checked by comparing the energy of the TP5 lattice 
previously described with that of "scarred" configurations obtained by decorating 
the original toroid in such a way that each +1 disclination on the external equator is 
replaced by a 5 — 7 — 5 mini-scar. The result of this comparison is summarized in the 
phase diagram of Fig. 53 in terms of r and the number of vertices of the triangular 
lattice V (the corresponding hexagonal lattice has twice the number of vertices, 
i.e. Vhex = 2y). V can be derived from the angular separation of neighboring 
disclinations in the same scar by approximating V ~ A/ Ay, with Ay = the 
area of a hexagonal Voronoi cell and o the lattice spacing. When the aspect ratio 
is increased from 1 to 6.8 the range of the curvature screening becomes shorter 
and the number of subunits required to destroy the stability of the TP5 lattice 
decreases. For r > 6.8, however, the geodesic distance between the two equators of 
the torus becomes too small and the repulsion between like-sign defects takes over. 
Thus the trend is inverted. 

Some example of toroidal lattices obtained from numerical simulations are shown 



in Fig. 54 The configurations are found by optimizing a system of V point-like 
particles constrained to lie on the surface of a torus and interacting via a pair 
potential of the form Uij = \j\xi — Xj^ where | • | denotes the Euclidean distance in 
R^. The choice of the cubic potential is motivated here by the so called "poppy seed 
bagel theorem" [80] , according to which the configuration of points that minimizes 
the Riesz energy E = ^i^j l/I^Kj — Xj\^ on a rectifiable manifold of Hausdorff 
dimension d is uniformly distributed on the manifold for s > d. In the case of a 
torus of revolution this implies that for small s the points are mostly distributed 
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Figure 53. (Color online) Isolated defects and scar phases in the (r, V) plane. When the number of vertices 

V increases the range of the screening curvature becomes smaller than one lattice spacing and disclinations 
appear delocalized in the form of a 5 — 7 — 5 grain boundary mini-scar. 

on the exterior of the torus (the interior becomes completely empty in the limit 

V — > oo). As s is increased, however, the points cover a progressively larger portion 
of the sm'face. The distribution becomes uniform for s > 2. On the other hand, since 
the number of local minima of the Riesz energy increases with s, it is practical to 
choose a value not much larger than two. The choice 5 = 3 has the further advantage 
of modelling a real physical system of neutral colloidal particles assembled at an 
interface [36 and is therefore suitable for direct comparison with experiments on 
colloidal suspensions. The lowest energies found, as well as the number of defects in 
the corresponding configuration, are reported in Table [7| The complete set of data 
produced in these simulations together with a collection of interactive 3D graphics 
for each low energy configuration studied can be found on-line jl96j . 

5.5.3. The Fat Torus Limit 

We have seen that disclination defects, forbidden in the lowest energy state of 
a planar crystal, may be energetically favored on a substrate of non-vanishing 
Gaussian curvature. It is therefore natural to ask whether large curvature can 
completely destroy crystalline order by driving the proliferation of a sufficiently 
high density of defects. The resulting state would be amorphous. The problem 
of generating amorphous structures by tiling a two-dimensional curved space with 
identical rigid subunits has drawn attention over the years, particularly through the 
connection to the structure of such disordered materials as supercooled liquids and 
metallic glasses. Since the work of Frank [197 the notion of geometrical frustration 
arises frequently in investigations of supercooled liquids and the glass transition. A 
paradigmatic example is represented by the icosahedral order in metallic liquids and 
glasses which, although locally favored, cannot propagate throughout all of three- 
dimensional Euclidean space. A two-dimensional analog, consisting of a liquid of 
monodisperse hard disks in a 2-manifold of constant negative Gaussian curvature 
(the hyperbolic plane) was first proposed by Nelson and coworkers in 1983 [T7 1I198] . 
In such a system the impossibility of covering the entire manifold with a 6-fold 
coordinated array of disks mimics many aspects of the geometrical frustration of 
icosahedral order in three dimensions. In all these models of geometrical frustration, 
however, the origin of the disorder is primarily due to the short-range nature of the 
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(3,32) (3,35) (3,60) (3,68) (3,126) 




(3,130) (3,180) (3,220) (3,340) (3,420) 




(4,117) (4,119) (4,120) (4,121) (4,125) 




(4,126) (4,260) (6,87) (6,108) (6,112) 




(6,116) (6,180) (6,220) (6,340) (6,420) 




(6,500) (6,1000) (20,180) (20,220) (20,300) 

Figure 54. (Color online) Selected low energy configurations for toroidal lattices of aspect ratio r = 3, 4, 6 
and 20. Lattices are labeled by (r, V), with V the number of particles. 



potential between the subunits. In a more realistic setting, part of the frustration 
is relieved by the fact that hexagonal unit cells can compress in order to match the 
underlying geometry. 

The embedding of a triangular lattice on an axisymmetric torus, provides a 
particularly suitable playground to study curvature-driven disorder. When r — > 1 
the Gaussian curvature on the inside of the torus grows like l/(r — 1) and diverges 
on the internal equator at ^/^ = vr. We thus expect a high density of defects in the 
vicinity of the curvature singularity and a resultant loss of the local 6— fold bond 
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Tabic 7. Low energy configuration for a selected number of 
toroidal lattices with aspect ratios r — 3, 4, 6 and 20. For each as- 
pect ratio the table displays the number of particles V, the lowest 
energy found and the number of fc— fold vertices Vk with k — 4-8. 
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orientational order. In this regime the system will have crystalline regions on the 
outside of the torus and amorphous regions near the curvature singularity. 

In this section we this claim can be supported using the elastic theory of contin- 
uous distributions of edge dislocations on a "fat" torus. Our argument is based on 
the following construction. As a consequence of the curvature singularity the sur- 
face area of an arbitrary wedge of angular width A(j) becomes smaller and smaller 
as the sectional angle ip increases and vanishes at V = ''J"- If a defect-free lattice 
is embedded on such a wedge, Bragg rows will become closer and closer as the 
singularity is approached with a consequent rise in the elastic energy (see Fig. 55). 
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Figure 55. (Color online) Top view of a defect free triangulation of a fat torus with (n, m, I) = (10, 10, 20) 
and V = 400. The corresponding elastic energy becomes very high in the interior of the torus where the 
triangles are more compressed to match the reduction of surface area. 



An intuitive way to reduce the distortion of the lattice is to recursively remove 
Bragg rows as one approaches the point tp = tt (see Fig. 56). This is equivalent 
to introducing a growing density of edge dislocations. This dislocation "cloud" 
will ultimately disorder the system by destroying the local 6— fold bond orienta- 
tional order. One might therefore view the curvature as playing the role of a local 
effective temperature which can drive "melting" by liberating disclinations and 
dislocations. In two-dimensional non-Euclidean crystals at T = 0, however, the 
mechanism for dislocation proliferation is fundamentally different from the usual 
thermal melting. While the latter is governed by an entropy gain due to unbinding 
of dislocation pairs, the amorphization at T = is due to the adjustment of the 
lattice to the geometry of the embedding manifold via the proliferation of defects 
and the consequent release of elastic stress. A similar phenomenon occurs in the 
disorder-driven amorphization of vortex lattices in type-II and high- Tc supercon- 
ductors [T99i[2nni[2nT] . 




Figure 56. (Color online) A schematic example of a dislocation pile-up on a square lattice resulting from 
the shrinking of the area on a regular wedge of a fat torus. 



Since the shrinking area per plaquette on the inside of the torus necessitates 
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a high density of dislocations we may approximate the dislocation cloud in this 
region by a continuous distribution of Burgers vector density b. Minimizing the 
elastic energy with respect to b yields a variational equation from which the optimal 
dislocation density can be calculated as a function of the ratio ed/iYE?) between 
the dislocation core energy and elastic energy scale YB?' with R = Ri = R2- 

As a starting point, we calculate the Green function GL{x,y) in the fat torus 
limit r ^ 1 (i.e. k ^ oo and u; — > 0, see Eq. (174])). The conformal angle ^ in this 
limit is 



Imi ^ = tan - , 



and to leading order of k we have 



^ / 2 



167r2 



— tan — 
47r 2 



^ Re{Li2(ae^'/')} ^ + ^ log (cos ^ 



To handle the limit of the Jacobi theta function we can take u = IS.z/k^ q = 

27r 

giTTT — g — g^j^^ calculate the limit q ^ 1. This can be done by using the modular 
transformation properties of Jacobi functions ^165| 1166] : 

^1 (^1 - ^) = -ii-iry^e^M^lr) ■ (218) 
Thus t' = — 1/t = iK/2, u' = u/t = Az/(2i) and q' = e"' = e~^, where 



lim ??i(ti, q) = lim « — e^^' 'di{u ,q) . 

This is easily evaluated by means of the expansion 

'&i{u, q) =2(74 sin u + o iq~^ 
Taking the logarithm and neglecting irrelevant constant terms, we obtain 



log 



^1 



z — z 



2i 

K 



log 



sinh 



z — z 



which finally leads to 

Gl(V,</>,V'','A') 



■4vr^*""T + 2^^°^ 



sinh 



z — z 
2 



(219) 



with z = tan('(/;/2) + With the Green function in hand, we can calculate the 
effect of the curvature singularity ai ^lJ = tt on the distribution of defects. Let b be 
the Burgers vector density of the dislocation cloud. Hereafter we work in a local 
frame, so that 



b = h^g^ + h^g^ , 



(220) 
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with Qi = diR a basis vector in the tangent plane of the torus, is a three- 
dimensional surface vector parameterizing the torus and its Cartesian components 



are given in Eq. (165). The quantity h has to be such that 



/ Sx h{x) = br, 
Jd 



with b£) the total Burger's vector in a generic domain D. Because on a closed man- 
ifold dislocation lines cannot terminate on the boundary, extending the integration 
to the whole torus we have: 



(221) 



Since the basis vectors Qi in Eq. (220) have the dimension of length, contravariant 



coordinates 6* have dimensions of an inverse area. Assuming all defects to be paired 
in the form of dislocations (i.e. qi = everywhere), the total energy of the crystal 
reads 



F = ^ y d^xT\x) + edj <fx , 



(222) 



where is the dislocation core energy. The function T{x) encoding the elastic 
stress due to the curvature and the screening contribution of the dislocation cloud 
obeys 



1 

Y 



^gT{x) = eiVib'^ix) - K{x) 



(223) 



where Vj is the usual covariant derivative along the coordinate-direction i and 
is the Levi-Civita antisymmetric tensor on the torus: 



gike 



.jk 



The stress function r(a;) can be expressed in the form T{x) = Td{^p,(f)) — Ts{iIj) 
with 



TsW 
Y 

^s(V^^. 

Y 



log 



1 



+ 1, 



2(1 + cosV')_ 
d^yeiVib\y)GL{x,y). 



(224a) 
(224b) 



The integral term in Eq. (224b) can now be reduced to a more friendly functional 



of b, suitable for a variational approach. Given the azimuthal symmetry we assume 
that all dislocations are aligned along b = b'^gs- 



Y 



1 

2^ 



d'^'^/gb'^i'ip') + sin -0' + vr sgn^ip — "0')] • 



(225) 



Substituting Eq. (225) and (224a) in Eq. (222) and minimizing with respect to b^, 
the problem can be converted, after some algebraic manipulations, to the following 
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Figure 57. (Color online) The Burgers vector component for different choices of A. 

Fredholm equation of the second kind: 



(226) 



where A = ed/iYR^), B{ij) = R^{l + cos^ljfh't'{il)) and the kernel JC{ip,ip') is given 
by: 



/C(V', V') = TT^ ^ < vr~i(^' + sin i;'){iP + sin V') 

4(1 + cosip ) I 

+ \'>p — ip \ + 2 cos — - — sm > . (227) 



The function /(V') on the right hand side of Eq. (226) is given by 

I J 1 1 



m 



dV'(i + cosVOr.(V'') 



= - { [log 2(1 + cos V') - 2] sin - 2 Ch{^p + tt) } (228) 
where CI2 is the Clausen function (see Ref. [202], pp. 1005-1006) defined as 

Chix) = -[dx log (2 sin = f; . 

As previously noted the dislocation core energy is is much smaller than the 
elastic energy scale YB?'. Eq (226) is then suitable to be solved in powers of the 
dimensionless number A: 

fi(V) = 5o(V') + \Bi{^) + \^B2{^) + ■■■ 

The corrections to the zero-order term Bq(iP) can be calculated recursively by 
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solving a set of Fredholm equations of the first kind: 



The function Bq(iP) associated with the Burgers vector density of the dislocation 
cloud in the limit A — > 0, on the other hand, can be calculated directly from 



Eq. (223) by setting the effective topological charge density on the right hand side 



to zero: 



eiVib\x)-Kix)=0. 
from which one obtains Bq{iP) = sinip and 



ii2(i + cosV')2 • 



(229) 



(230) 



The Burgers vector density b'^ obtained by a numerical solution of Eq. (226) is 



shown in Fig. 57 for different values of A. The Burgers vector density is measured 
in units of R~'^. The function b'^ has cubic singularities at ip = ibvr and is approxi- 



mately zero on the outside of the torus. The solid blue curve in Fig. 57 represents 



the zeroth order solution of Eq.(230). 



Now, in the theory of dislocation mediated melting a system at the solid liquid 
phase boundary is described as a crystalline solid saturated with dislocations. In 
three-dimensions, in particular, there is a strong experimental evidence of the exis- 
tence of a critical dislocation density at the melting point p{Tm) ~ 0.66"^ where b 
is the length of the length of the smallest perfect-dislocation Burgers vector |203] . 
Several theoretical works have motivated this evidence both for three-dimensional 
solids and vortex lattices in super conductors [200] . On the other hand, given the 
existence of such a critical density, its value can be empirically used to determine 
whether a system is in a solid or liquid-like phase in the same spirit as the Lin- 
demann criterion. With this goal in mind we can calculate the dislocation density 
by requiring |&| = pa with p the density of single lattice spacing dislocations. This 
yields: 



V3 



V 



tan 



+ o(A) 



(231) 



58 



Solving p{ip, V)a'^ = 0.6 as a function of -0 and V we obtain the diagram of Fig. 
As expected the inside of the torus contains an amorphous region whose angular 
size decreases with the number of vertices F as a consequence of the reduction of 
the lattice spacing. 



6. Conclusion and discussion 

In this article we have reviewed recent progress toward understanding the ground 
state properties of two-dimensional ordered phases on substrate of non-zero Gaus- 
sian curvature. Curved phases of matter are found in a broad class of systems com- 
prising materials as different as viral capsids and carbon nanotori. The geometry 
of the substrate gives rise to novel defective structures that would be energetically 
prohibitive in the ground state of conventional flat systems. Defects may appear 
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Figure 58. (Color online) Phase diagram for curvature driven amorphization. The inside of the torus 
contains an amorphous region whose angular size decreases with the number of vertices y as a consequence 
of the reduction of the lattice spacing. 




Figure 59. (Color online) (Left) Idealized drawing of a side and top view and of a rippled gold nanoparticle 
showing the two polar defects that must exist to allow the alternation of concentric rings. (Right) TEM 
images of a chain obtained when 11-mercaptoundecanoic acid functionalized nanoparticles are reacted with 
a water phase containing divalent 1,6-diaminohexane in a two-phase reaction. Scale bar 50 nm. From |204| . 

in curved space either because they are required by the topology of the underlying 
substrate (as in the case of spherical crystals and nematics) or because they are 
favored by the curvature itself. The purely curvature driven energetic nature of the 
latter phenomenon is of fundamental interest. 

On the material science side, non-Euclidean systems provides a promising route 
to constructing arrays of nanoparticles via the chemical functionalization of topo- 
logical defects created on the surface of the particles by coating them with an or- 
dered monolayer. Candidates for such coatings include triblock coplymers, gemini 
lipids, metallic or semiconducting nanorods and conventional liquid crystal com- 
pounds. If the induced order is vector-like or striped, then there must necessarily 
be two defective regions corresponding to the bald spots exhibited by a combed 
sphere or the source and sink of fluid flow confined to a sphere. Not only are the 
defects physically and mathematically distinguished, since the condensed matter 
order vanishes there, but it also turns out that chemical processes can detect the 
defects and insert linker molecules at the precise defect locations. This creates the 
possibility of efficiently making functionalized nanoparticles with a precise valence 
and corresponding directional bonding. The nanoparticles themselves can then be 
linked into well-defined two and three dimensional arrays. Control over the valence 
and the geometry of the directional bonding can be achieved by varying the nature 
of the ordered monolayer. 

This scheme has been demonstrated recently in the work of DeVries et al. [2D4J in 
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which gold nanoparticles are coated with two species of naturahy phase separating 
hgands (see Fig. 59). Phase separation on the spherical gold surface translates to a 
striped arrangement of alternating ligands. Such a spherical smectic has topological 
defects at the north and south pole of the nanopartcle. These may be functionalized 
by attaching thiol-based linker molecules at the defects via place-exchange reac- 
tions. This creates a divalent gold nanoparticle with directional bonds 180 degrees 
apart. These in turn can be linked to create polymers and free-standing films. 
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